R/L across languages: Modeling

1 Introduction

This script uses the output of r_l_preparation.Rmd and produes all the values and plots reported in the paper.

if( !all(file.exists("./data/web_raw_trials.csv", "./data/web_experiment_cleaned.csv", "./data/field_raw_trials.csv", "./data/field_experiment_cleaned.csv", "./2_r_l_preparation.html")) )
{
  # Seems like the script 2_r_l_preparation.Rmd was not invoked, so let's compile it!
  require(rmarkdown);
  #rmarkdown::render("./2_r_l_preparation.Rmd");
  xfun::Rscript_call(rmarkdown::render, list(input="./2_r_l_preparation.Rmd"));
  if( !all(file.exists("./data/web_raw_trials.csv", "./data/web_experiment_cleaned.csv", "./data/field_raw_trials.csv", "./data/field_experiment_cleaned.csv", "./2_r_l_preparation.html")) )
  {
    # Seems this issue is quite serious?
    stop("Failed to compile the '2_r_l_preparation.Rmd' Rmarkdonw script: please try to do manually!\n");
  }
}

2 Setup

# Load packages:

library(tidyverse) # data processing
library(brms) # bayesian models
#library(cmdstanr) # install it with: install.packages("cmdstanr", repos = c("https://mc-stan.org/r-packages/", getOption("repos"))) followed by install_cmdstan()
library(ggdist) # for plotting
library(tidybayes) 

# option for Bayesian regression models:
# use all available cores for parallel computing
options(mc.cores = parallel::detectCores())

## Set the script's path as working directory
#parentfolder = rstudioapi::getActiveDocumentContext()$path 
#setwd(dirname(parentfolder))
#parentfolder <- getwd()
parentfolder <- "."; # assume current folder is the document folder

#models        <- paste0(parentfolder, '/models/')
plots         <- paste0(parentfolder, '/plots/')
data          <- paste0(parentfolder, '/data/')

if( !dir.exists("./cached_results") ) dir.create("./cached_results", showWarnings=FALSE);

# Various auxiliary functions:
library(parallel);
library(lme4);
library(performance);
library(brms);
library(bayestestR);
library(ggplot2);
library(gridExtra);
library(ggpubr);
library(sjPlot);
library(knitr);

brms_ncores  <- max(detectCores(all.tests=TRUE, logical=FALSE), 4, na.rm=TRUE); # try to use multiple cores, if present

# Log odds to probability (logistic regression intercept):
lo2p <- function(x){ o <- exp(x); return (o/(1+o));}

# Log odds to odds ratio change as percent (logistic regression slopes):
lo2or <- function(x){ o <- exp(x); if(x < 0){ return (-(1-o)) } else { return (o-1) };}

# Log odds to odds ratio change between level values (logistic regression slopes):
#lo2ps <- function(a,b){if(b > 0){ return (plogis(a+b)-plogis(a)) }else{ return (-(plogis(a)-plogis(a+b))) }}
lo2ps <- function(a,b){ plogis(a+b) - plogis(a) }

# Scientific notation using Markdown conventions (inspired from https://www.r-bloggers.com/2015/03/scientific-notation-for-rlatex/):
scinot <- function(xs, digits=2, pvalue=TRUE)
{
  scinot1 <- function(x)
  {
    sign <- "";
    if(x < 0)
    {
      sign <- "-";
      x <- -x;
    }
    exponent <- floor(log10(x));
    if(exponent && pvalue && exponent < -3) 
    {
      xx <- round(x / 10^exponent, digits=digits);
      e <- paste0("×10^", round(exponent,0), "^");
    } else 
    {
      xx <- round(x, digits=digits+1);
      e <- "";
    }
    paste0(sign, xx, e);
  }
  vapply(xs, scinot1, character(1));
}

# Escape * in a string:
escstar <- function(s)
{
  gsub("*", "\\*", s, fixed=TRUE);
}

# Figure and Table caption adapted from https://stackoverflow.com/questions/37116632/rmarkdown-html-number-figures: 
outputFormat = opts_knit$get("rmarkdown.pandoc.to"); # determine the output format of the document
if( is.null(outputFormat) ) outputFormat = ""; # probably not run within knittr
capTabNo = 1; capFigNo = 1; # figure and table caption numbering, for HTML do it manually
#Function to add the Table Number
capTab = function(x)
{
  if(outputFormat == 'html'){
    x = paste0("**Table ",capTabNo,".** ",x,"")
    capTabNo <<- capTabNo + 1
  }; x
}
#Function to add the Figure Number
capFig = function(x, show_R_version=TRUE, show_package_versions=NULL, is_map=FALSE)
{
  if(outputFormat == 'html')
  {
    x <- paste0("**Figure ",capFigNo,".** ",x,"");
    if( show_R_version || (!is.null(show_package_versions) && length(show_package_versions) > 0) )
    {
      x <- paste0(x, " Figure generated using ");
      if( show_R_version ) x <- paste0(x, stringr::str_replace(R.version.string, stringr::fixed("R "), "[`R`](https://www.r-project.org/) "));
      if( !is.null(show_package_versions) && length(show_package_versions) > 0 )
      {
        x <- paste0(x, ifelse( show_R_version, " and ", " "));
        x <- paste0(x, ifelse( length(show_package_versions) > 1, "packages ", "package "));
        x <- paste0(x, paste0(vapply(show_package_versions, function(x) paste0("`",x,"`", " (version ", packageVersion(x),")"), character(1)), collapse=", "), ".");
      }
      if( is_map ) x <- paste0(x, " Maps are using public domain data from the [Natural Earth project](https://www.naturalearthdata.com/) as provided by the `R` package `maps`.");
    }
    capFigNo <<- capFigNo + 1;
  }; 
  x;
}

# For reproducible reporting (see also the endof this document):

packageVersion('tidyverse')
packageVersion('ggplot2')
packageVersion('brms')
#packageVersion('cmdstanr')
packageVersion('ggdist')

# Load ggplot2 theme and colors:

source('theme_timo.R')

colorBlindBlack8  <- c("#000000", "#E69F00", "#56B4E9", "#009E73", 
                       "#F0E442", "#0072B2", "#D55E00", "#CC79A7")

# Load data:

web     <- read_csv(paste0(data, 'web_experiment_cleaned.csv'))
web_raw <- read_csv(paste0(data, '/web_raw_trials.csv'))

field     <- read_csv(paste0(data, 'field_experiment_cleaned.csv'))
field_raw <- read_csv(paste0(data, '/field_raw_trials.csv'))

3 Online experiment

3.1 Descriptive statistics

First, how many participants?

nrow(web)

## [1] 903

Sex division

table(web$Sex)

## 
##   F   M 
## 681 222

Ages

summary(web$Age)

##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##   18.00   23.00   29.00   32.92   40.00   84.00

Order division

# Counts
table(web$Order)

## 
## l_first r_first 
##     436     467

# Percentage
prop.table(table(web$Order)) * 100

## 
## l_first r_first 
## 48.2835 51.7165

First, how many languages?

web %>% count(Name) %>% nrow()

## [1] 25

Does this number correspond with the L1s?

web %>% count(Language) %>% nrow()

## [1] 25

How many families?

web %>% count(Family) %>% nrow()

## [1] 9

How many have the R/L distinction in the L1 among the languages?

web %>% count(Language, r_l_distinction_L1) %>% count(r_l_distinction_L1)

How many really use the alveolar trill in L1 among the languages?

web %>% count(Language, trill_real_L1) %>% count(trill_real_L1)

How many really have the alveolar trill in L1 as an allophone among the languages?

web %>% count(Language, trill_occ_L1) %>% count(trill_occ_L1)

What about the same questions for L2. But this will not neatly sum up to 25, due to various possible scenarios for L2 within a specific L1.

How many have the R/L distinction in the L2 among the languages?

web %>% count(Language, r_l_distinction_L2) %>% count(r_l_distinction_L2)

How many really use the alveolar trill in L2 among the languages?

web %>% count(Language, trill_real_L2) %>% count(trill_real_L2)

How many really have the alveolar trill in L2 as an allophone among the languages?

web %>% count(Language, trill_occ_L2) %>% count(trill_occ_L2)

What is the grand average congruent behavior?

mean(web$Match)

## [1] 0.8726467

87.3%

What about only among those who have L1 without the distinction?

web %>%
  filter(r_l_distinction_L1 == "0") %>%
  summarize(mean_match = mean(Match, na.rm = TRUE))

83.9%

What about only among those who have L1 without the distinction and no L2 that distinguishes?

web %>%
  filter(r_l_distinction_L1 == "0") %>%
  filter(!EnglishL2YesNo) %>% 
  filter(r_l_distinction_L1 == '0') %>% 
  summarize(mean_match = mean(Match, na.rm = TRUE))

85.1%

Compute average matching behavior per language and sort:

web_avg <- web %>%
  group_by(Language) %>% 
  summarize(M = mean(Match)) %>% 
  arrange(desc(M)) %>% 
  mutate(percent = round(M, 2) * 100,
         percent = str_c(percent, '%'))

# Show:

web_avg %>% print(n = Inf)

## # A tibble: 25 × 3
##    Language     M percent
##    <chr>    <dbl> <chr>  
##  1 FI       1     100%   
##  2 FR       0.982 98%    
##  3 EN       0.974 97%    
##  4 DE       0.953 95%    
##  5 SE       0.952 95%    
##  6 DK       0.944 94%    
##  7 HU       0.941 94%    
##  8 JP       0.927 93%    
##  9 KR       0.909 91%    
## 10 GR       0.905 90%    
## 11 PL       0.887 89%    
## 12 RU       0.872 87%    
## 13 EE       0.860 86%    
## 14 FA       0.857 86%    
## 15 AM       0.85  85%    
## 16 ZU       0.85  85%    
## 17 IT       0.846 85%    
## 18 TR       0.811 81%    
## 19 ES       0.806 81%    
## 20 GE       0.8   80%    
## 21 TH       0.8   80%    
## 22 PT       0.770 77%    
## 23 RO       0.742 74%    
## 24 AL       0.7   70%    
## 25 CN       0.696 70%

Check some demographics, also to report in the paper. First, the number of participants per language:

web %>% 
  count(Name, sort = TRUE) %>% print(n = Inf)

## # A tibble: 25 × 2
##    Name           n
##    <chr>      <int>
##  1 German        85
##  2 Portuguese    61
##  3 French        57
##  4 Japanese      55
##  5 Polish        53
##  6 Italian       52
##  7 Russian       47
##  8 Chinese       46
##  9 Estonian      43
## 10 Greek         42
## 11 English       39
## 12 Turkish       37
## 13 Spanish       36
## 14 Hungarian     34
## 15 Romanian      31
## 16 Korean        22
## 17 Farsi         21
## 18 Swedish       21
## 19 Armenian      20
## 20 Thai          20
## 21 Zulu          20
## 22 Danish        18
## 23 Finnish       18
## 24 Georgian      15
## 25 Albanian      10

Then, the number of L1 speakers who have R/L distinction vs. who don’t:

web %>% count(r_l_distinction_L1) %>%
  mutate(prop = n / sum(n),
         percent = round(prop, 2) * 100,
         percent = str_c(percent, '%'))

How many people do not have any L2?

sum(is.na(web$L2)) / nrow(web)

## [1] 0.1351052

# bilinguals:
1 - sum(is.na(web$L2)) / nrow(web)

## [1] 0.8648948

Check how many people knew English as their L2:

web %>% count(EnglishL2YesNo) %>% 
  mutate(prop = n / sum(n),
         percent = round(prop, 2) * 100,
         percent = str_c(percent, '%'))

Of those that don’t use a R/L distinction in their L1, how many use R/L distinction in their L2? (double-check if logic alright!)

web %>%
  filter(r_l_distinction_L1 == 0) %>% 
  count(r_l_distinction_L2 == 1) %>% 
  mutate(prop = n / sum(n),
         percent = round(prop, 2) * 100,
         percent = str_c(percent, '%'))

How many “pure” speakers were there?, i.e., those people that 1) don’t know English, 2) don’t use an L1 with a R/L distinction, and 3) don’t know an L2 that distinguishes R/L.

web %>% 
  filter(r_l_distinction_L1 == 0) %>%  # excludes English as well
  filter(!EnglishL2YesNo) %>% 
  filter(r_l_distinction_L2 == 0) %>% 
  nrow()

## [1] 1

1 person!

Let’s check if this is correct. This gives the list of all participants for whom this applies.

web %>% 
    filter(r_l_distinction_L1 == 0 & !EnglishL2YesNo & r_l_distinction_L2 == 0) %>% 
    print(n = 50)

## # A tibble: 1 × 18
##   ID         Match Language Sex     Age Name    Script Family Autotyp_Area L2   
##   <chr>      <dbl> <chr>    <chr> <dbl> <chr>   <chr>  <chr>  <chr>        <chr>
## 1 JP_2040343     1 JP       F        48 Japane… diffe… Japan… N Coast Asia chin…
## # ℹ 8 more variables: EnglishL2YesNo <lgl>, Order <chr>,
## #   r_l_distinction_L1 <dbl>, trill_real_L1 <dbl>, trill_occ_L1 <dbl>,
## #   r_l_distinction_L2 <dbl>, trill_real_L2 <dbl>, trill_occ_L2 <dbl>

Are these really “pure”? What languages do they speak?

web %>% 
  filter(r_l_distinction_L1 == 0) %>%  # excludes English as well
  filter(!EnglishL2YesNo) %>% 
  filter(r_l_distinction_L2 == 0) %>% 
  count(Language)

One Japanese speaker.

Nevertheless, let’s explore whether these also show matches?

web %>% 
  filter(r_l_distinction_L1 == 0) %>%  # excludes English as well
  filter(!EnglishL2YesNo) %>% 
  filter(r_l_distinction_L2 == 0) %>% 
  count(Match) %>% 
  mutate(prop = n / sum(n),
         percent = round(prop, 2) * 100,
         percent = str_c(percent, '%'))

Yes, similar to above 85%.

Now, some relevant plots:

Distribution of matches in the web data by Language.

Distribution of matches in the web data by Family.

Distribution of matches in the web data by Area.

Distribution of matches in the web data by [r]/[l] distinction in L1.

Distribution of matches in the web data by [r] is main allophone in L1.

Distribution of matches in the web data by [r]/[l] distinction in L2.

Distribution of matches in the web data by [r] is main allophone in L2.

Distribution of matches in the web data (i.e., not modelled but the actual raw data) by Order (columns) and the [r] is main allophone in the L1 (rows).

3.2 What is the chance level?

Our experiment has two trials and counterbalancing of order, but we code only strict matches, meaning that there are four combination of responses for two values of Match: 1/1 = match, 1/0, 0/1 = partial matches (i.e., no match), and 0/0 = mismatch (i.e., no match). If we assume that that the choices for the two sounds are independent, then We can simulate this:

sim_chance_level <- replicate(1000, # simulate 1000 "experiments"
                              mean(replicate(1000, # the proportion of "matches"  for 1000 random "participants"
                                             ("r" == sample(c("r","l"), size=1, prob=c(0.5, 0.5))) && 
                                               ("l" == sample(c("r","l"), size=1, prob=c(0.5, 0.5)))))); # only full matches matter (here order does not matter as the repose is random
hist(100*sim_chance_level, main="Histogram of chance level", xlab="% full matches", col="gray90", xlim=c(0,100));
abline(v=100*mean(sim_chance_level), col="blue", lwd=2);

Figure 1. Simulating the % of ‘full’ matches across 1000 ‘experiments’ each with 1000 ‘participants’ that associate randomly (at 50%) [r] and [l] with the wavy or the straight line, assuming independence between the two choices. Figure generated using R version 4.3.3 (2024-02-29)

Indeed, as expected, the chance level of such a ‘full’, conservative match is indeed 25%.

However, it is clear that we cannot assume independence, as the choice made for the first sound clearly affects the choice made for the second, but it is unclear how much.

At the other extreme, if we assume perfect dependence between the two choices (i.e., whatever is chosen for the first sound, will not be chosen for the second), we have the following baseline probability:

sim_chance_level2 <- replicate(1000, # simulate 1000 "experiments"
                              mean(replicate(1000, # the proportion of "matches"  for 1000 random "participants"
                                             ("r" == sample(c("r","l"), size=1, prob=c(0.5, 0.5)))))); # only the first choice really matters (the second is fully dependent on it, so if the first is correct, the second will be as well, resulting in a match, and vice-versa)
hist(100*sim_chance_level2, main="Histogram of chance level", xlab="% full matches", col="gray90", xlim=c(0,100));
abline(v=100*mean(sim_chance_level2), col="blue", lwd=2);

Figure 2. Simulating the % of ‘full’ matches across 1000 ‘experiments’ each with 1000 ‘participants’ that associate randomly (at 50%) [r] and [l] with the wavy or the straight line, assuming perfect dependence between the two choices. Figure generated using R version 4.3.3 (2024-02-29)

Indeed, as expected, the chance level of such a ‘full’, conservative match is indeed 50%.

It is clear that we are in between these extremes, probably closer to 50% than to 25%, so we will use the 50% chance level throughout being very clear that this is a conservative level.

([)Please note that we abstain from attempting to estimating the chance level from the data, as it is very hard to disentangle, in our experimental design, the chance level from the sound symbolism bias we want to estimate.)

3.3 Regression models

3.3.1 Making sense of the fixed effects

In the various logistic regression models that we fit here, we use various techniques for judging the contribution of any given fixed effect:

• point estimates and 95% Credible Intervals (CIs; see bayestestR::ci with method = "ETI" for details), which should not, in principle, include 0,
• plots of the posterior distribution of the parameters; these should also not include 0,
• formal hypothesis tests against 0 (see brms::hypothesis for details): these are always directional and in the same direction as the point estimate, and should not be taken, except in the case of the intercept, as being a priori justified, but simply as a way of summarizing where the bulk of the posterior distribution is relative to 0; we use the default cut-off of 0.95 to judge the “existence” of the effect in the given direction (marked with a star “*“).

3.3.2 What predictors to include?

Here we are interested in predicting the probability of a perfect match, Match, from a set of potential predictors and their interactions. A priori, given our hypothesis, the potential predictors here are the order of presentation, Order, and the various properties of the L1(s) and L2(s) spoken by the participant in what concerns [r] and [l], namely:

• if [r] and [l] are phonologically distinct in any of the L1(s), r_l_distinction_L1, or L2(s), r_l_distinction_L2, respectively,
• if [r] is the main allophone in any of the L1(s), trill_real_L1, or L2(s), trill_real_L2, respectively, and
• if [r] appears as an allophone at all in some lect or register of the L1(s), trill_occ_L1, or L2(s), trill_occ_L2, respectively.

Importantly, we have no a priori reason to consider the participant’s Sex or Age as potential predictors, but we can nevertheless test their influence on Match.

We code here all the [r] and [r]/[l] predictors as factors with a treatment contrast as follows:

• r_l_distinction_L1 and r_l_distinction_L2 can be “same” (baseline) or “distinct”,
• trill_real_L1, trill_real_L2, trill_occ_L1 and trill_occ_L2 can be “no” (baseline) or “yes”.

On the other hand, we code Order as a factor with possible values “r_first” and “l_first” but using a contrast (or deviation) coding (-0.5, +0.5), because Order is (almost) balanced at 51.7% “r_first” vs 48.3% “l_first”, which ensures that the intercept represents (roughly) the grand mean (see, for example, Chapter Sum contrasts in An Introduction to Bayesian Data Analysis for Cognitive Science by Bruno Nicenboim, Daniel Schad, and Shravan Vasishth). (Please note that the other predictors are not balanced, making the treatment contrast better suited.)

# Make sure we code the factors with the intended baseline levels and contrasts:
web$Order              <- factor(web$Order, levels = c("r_first", "l_first")); contrasts(web$Order) <- c(-0.5, +0.5)
web$r_l_distinction_L1 <- factor(c("same", "distinct")[web$r_l_distinction_L1 + 1], levels=c("same", "distinct"));
web$trill_real_L1      <- factor(c("no", "yes")[web$trill_real_L1 + 1], levels=c("no", "yes"));
web$trill_occ_L1       <- factor(c("no", "yes")[web$trill_occ_L1 + 1], levels=c("no", "yes"));
web$r_l_distinction_L2 <- factor(c("same", "distinct")[web$r_l_distinction_L2 + 1], levels=c("same", "distinct"));
web$trill_real_L2      <- factor(c("no", "yes")[web$trill_real_L2 + 1], levels=c("no", "yes"));
web$trill_occ_L2       <- factor(c("no", "yes")[web$trill_occ_L2 + 1], levels=c("no", "yes"));

web %>% count(Match) %>% mutate(prop = paste0(round(100* n / sum(n),1),"%")) 

# match is quite unbalanced: 12.7% vs 87.3%

web %>% count(Order) %>% mutate(prop = paste0(round(100* n / sum(n),1),"%")) 

# the order effect is decently balanced: 51.6% vs 48.3%

web %>% count(r_l_distinction_L1) %>% mutate(prop = paste0(round(100* n / sum(n),1),"%")) 

# highly imbalanced: 15.8% vs 84.2%

web %>% count(trill_real_L1) %>% mutate(prop = paste0(round(100* n / sum(n),1),"%")) 

# ok: 58.8% vs 41.2%

web %>% count(trill_occ_L1) %>% mutate(prop = paste0(round(100* n / sum(n),1),"%")) 

# 100% are 1 --> excluded from the model

## And for L2, just in case
web %>% count(r_l_distinction_L2) %>% mutate(prop = paste0(round(100* n / sum(n),1),"%")) 

# 13.5% missing, the rest almost all (86.3%) are 1 and 0.2% are 0 ---> exclude as well

web %>% count(trill_real_L2) %>% mutate(prop = paste0(round(100* n / sum(n),1),"%")) 

# 13.5% missing, the rest is balanced: 53.8% vs 32.7%

web %>% count(trill_occ_L2) %>% mutate(prop = paste0(round(100* n / sum(n),1),"%")) 

# 13.5% missing, the rest are all (86.5%) are 1 ---> exclude as well

It can be seen that, in our data, trill_occ_L1 and trill_occ_L2 are both 100% “yes” (for the non-missing data cases), which means that we will not include them in any further analyses.

3.3.3 Note about logistic regression coefficients

Because we want to predict the probability of a Match, which is a binary variable, we will use logistic regression throughout. However, logistic regression coefficients are notoriously hard to interpret, as they are reported on the log odds scale. Here, we show the intercepts and regression slopes both as log odds ratios and as probabilities, as appropriate. (See, for example, here, here and here for more explanations.)

Please note that for the intercept (α) this represents the probability of a match when all the predictors are at 0.0 or their baseline levels; for example, an intercept of 2.87 corresponds to a probability of 94.6% of a match when holding all the other predictors at 0 or at base level.

For the slopes (β), this interpretation changes, and we show both the equivalent percent change in odds and the change in the probability of a match between this predictor’s baseline and non-baseline levels (or for a unit increase from 0.0 to 1.0) when all the other predictors are at baseline (or at 0.0). For example, a slope of -0.84 for a binary DV corresponds to an odds ratio of 0.43, which represents a decrease of 56.8% in the odds of a match, or, equivalently, a decrease by 6.2% in the probability of a match from the baseline 94.6% when the DV is not at its baseline level to 88.4%.

3.3.4 What random structure to use?

An important point to clarify is what should be the random effects structure of our regression models.

A priori, there are three potential “grouping factors”: Language, Family and (AUTOTYP) Area, that are meaningfully considered as random effects. We have 25 L1 languages, from 9 families and 6 areas:

AUTOTYP areas.* Reproduced from The AUTOTYP database GitHub repo under a CC-BY-4.0 license.

unique(web[,c("Name", "Family", "Autotyp_Area")]) %>% arrange(Autotyp_Area, Family, Name)

However, when it comes to the r/l distinction and the [r] is main allophone in the language (L1 or L2), it is clear that these variables do not, by definition, vary within Language (so there should be no random slope here), and, at least in our current data, they very very little within the levels of the other two grouping factors as well (see below), raising the legitimate question whether we should model random slopes for these two variables at all.

Order: we know this varies within all levels of Language, Family and Area because of the experimental design, so Order should have varying slopes for all three.

r/l distinction in L1 (r_l_distinction_L1):

table(web$r_l_distinction_L1, web$Language)

##           
##            AL AM CN DE DK EE EN ES FA FI FR GE GR HU IT JP KR PL PT RO RU SE TH
##   same      0  0 46  0  0  0  0  0  0  0  0  0  0  0  0 55 22  0  0  0  0  0  0
##   distinct 10 20  0 85 18 43 39 36 21 18 57 15 42 34 52  0  0 53 61 31 47 21 20
##           
##            TR ZU
##   same      0 20
##   distinct 37  0

table(web$r_l_distinction_L1, web$Family)

##           
##            Benue-Congo Finno-Ugric  IE Japanese Kartvelian Korean Sino-Tibetan
##   same              20           0   0       55          0     22           46
##   distinct           0          95 593        0         15      0            0
##           
##            Tai-Kadai Turkic
##   same             0      0
##   distinct        20     37

table(web$r_l_distinction_L1, web$Autotyp_Area)

##           
##            Europe Greater Mesopotamia Inner Asia N Coast Asia S Africa
##   same          0                   0          0           77       20
##   distinct    539                  93        108            0        0
##           
##            Southeast Asia
##   same                 46
##   distinct             20

This is perfectly uniform within the levels of each of the three factors, so random slopes are arguably not justified.

r/l distinction in L2 (r_l_distinction_L2):

table(web$r_l_distinction_L2, web$Language)

##           
##            AL AM CN DE DK EE EN ES FA FI FR GE GR HU IT JP KR PL PT RO RU SE TH
##   same      0  0  0  0  0  0  1  0  0  0  0  0  0  0  0  1  0  0  0  0  0  0  0
##   distinct  9 20 30 84 16 43 28 34 18 18 46 15 42 26 50 28 21 52 42 30 43 19 18
##           
##            TR ZU
##   same      0  0
##   distinct 29 18

table(web$r_l_distinction_L2, web$Family)

##           
##            Benue-Congo Finno-Ugric  IE Japanese Kartvelian Korean Sino-Tibetan
##   same               0           0   1        1          0      0            0
##   distinct          18          87 533       28         15     21           30
##           
##            Tai-Kadai Turkic
##   same             0      0
##   distinct        18     29

table(web$r_l_distinction_L2, web$Autotyp_Area)

##           
##            Europe Greater Mesopotamia Inner Asia N Coast Asia S Africa
##   same          1                   0          0            1        0
##   distinct    478                  82        104           49       18
##           
##            Southeast Asia
##   same                  0
##   distinct             48

There is almost no “same” at all, so random slopes are arguably not justified.

[r] is main allophone in L1 (trill_real_L1):

table(web$trill_real_L1, web$Language)

##      
##       AL AM CN DE DK EE EN ES FA FI FR GE GR HU IT JP KR PL PT RO RU SE TH TR
##   no   0  0 46 82 18  0 38  0  0  0 56 14 42  0  0 55 22  0 61  0  0 20 20 37
##   yes 10 20  0  3  0 43  1 36 21 18  1  1  0 34 52  0  0 53  0 31 47  1  0  0
##      
##       ZU
##   no  20
##   yes  0

table(web$trill_real_L1, web$Family)

##      
##       Benue-Congo Finno-Ugric  IE Japanese Kartvelian Korean Sino-Tibetan
##   no           20           0 317       55         14     22           46
##   yes           0          95 276        0          1      0            0
##      
##       Tai-Kadai Turkic
##   no         20     37
##   yes         0      0

table(web$trill_real_L1, web$Autotyp_Area)

##      
##       Europe Greater Mesopotamia Inner Asia N Coast Asia S Africa
##   no     317                  51          0           77       20
##   yes    222                  42        108            0        0
##      
##       Southeast Asia
##   no              66
##   yes              0

This is almost uniform within Language, but there is variation for the IE level of Family (almost 50%:50% between “no” and “yes”), and between 2 levels of Area (Europe with 60%:40% and Greater Mesopotamia with 55%:45% “no”:“yes”), but this is not enough to justify random slopes either.

[r] is main allophone in L2 (trill_real_L2):

table(web$trill_real_L2, web$Language)

##      
##       AL AM CN DE DK EE EN ES FA FI FR GE GR HU IT JP KR PL PT RO RU SE TH TR
##   no   5  1 29 38 13  5 16 22 16  8 28  3 28 18 30 26 21 34 19 21 32 11 18 28
##   yes  4 19  1 46  3 38 13 12  2 10 18 12 14  8 20  3  0 18 23  9 11  8  0  1
##      
##       ZU
##   no  16
##   yes  2

table(web$trill_real_L2, web$Family)

##      
##       Benue-Congo Finno-Ugric  IE Japanese Kartvelian Korean Sino-Tibetan
##   no           16          31 314       26          3     21           29
##   yes           2          56 220        3         12      0            1
##      
##       Tai-Kadai Turkic
##   no         18     28
##   yes         0      1

table(web$trill_real_L2, web$Autotyp_Area)

##      
##       Europe Greater Mesopotamia Inner Asia N Coast Asia S Africa
##   no     283                  48         45           47       16
##   yes    196                  34         59            3        2
##      
##       Southeast Asia
##   no              47
##   yes              1

Here there is enough variation within the levels for all three factors, justifying the inclusion of random slopes.

Given these, we did the following:

• we always include random slopes for Order by Language, Family and Area,
• we do not include random slopes for the r/l distinction in L1 (r_l_distinction_L1) or L2 (r_l_distinction_L2), nor for the presence of [r] in L1 (trill_real_L1), but
• we do include random slopes for presence of [r] in L2 (trill_real_L2) by Language, Family and Area.

3.3.5 Model 1: what is the overall probability of a match?

if( !file.exists("./cached_results/b1_res_web.RData") )
{
  # prior predictive checks:
  # what priors we can set (use Order for that):
  get_prior(Match ~ 1 + Order +
              (1 + Order | Language) +
              (1 + Order | Family) +
              (1 + Order | Autotyp_Area),
            data=web,
            family=bernoulli(link='logit'));
  # -> "sd" priors seem alright (student_t(3, 0, 2.5)), as does the "Intercept" (student_t(3, 0, 2.5)), but lkj(1) for "cor" might too accepting of extreme correlations, and (flat) for "b" is clearly not ok
  # so, we keep "sd" and "Intercept" but use lkj(2) for "cor" and student_t(5, 0, 2.5) for "b"
  b_priors <- brms::brm(Match ~ 1 + Order +
                          (1 + Order | Language) +
                          (1 + Order | Family) +
                          (1 + Order | Autotyp_Area),
                        data=web,
                        family=bernoulli(link='logit'),
                        prior=c(brms::set_prior("student_t(5, 0, 2.5)", class="b"), 
                                brms::set_prior("lkj(2)", class="cor")),
                        sample_prior='only',  # needed for prior predictive checks
                        seed=998, cores=brms_ncores, iter=10000, warmup=4000, thin=2, control=list(adapt_delta=0.999, max_treedepth=13));
  grid.arrange(pp_check(b_priors, type='stat', sta ='min') + xlab('p(match)') + ggtitle('Prior predictive distribution of minimum values'),
               pp_check(b_priors, type='stat', stat='mean') + xlab('p(match)') + ggtitle('Prior predictive distribution of means'),
               pp_check(b_priors, type='stat', stat='max') + xlab('p(match)') + ggtitle('Prior predictive distribution of maximum values'),
               ncol=1); # seems fine but our value is a bit extreme
  
  # the model:
  b1 <- brm(Match ~ 1 + Order +
              (1 + Order | Language) +
              (1 + Order | Family) +
              (1 + Order | Autotyp_Area),
            data=web,
            family=bernoulli(link='logit'),
            prior=c(brms::set_prior("student_t(5, 0, 2.5)", class="b"), 
                    brms::set_prior("lkj(2)", class="cor")),
            save_pars=save_pars(all=TRUE), # needed for Bayes factors
            sample_prior=TRUE,  # needed for hypotheses tests
            seed=998, cores=brms_ncores, iter=10000, warmup=4000, thin=2, control=list(adapt_delta=0.999, max_treedepth=13));
  summary(b1); mcmc_plot(b1, type="trace"); mcmc_plot(b1); # very decent
  mcmc_plot(b1, type="trace");
  bayestestR::hdi(b1, ci=0.95);
  hypothesis(b1, "Order1 = 0"); 
  # posterior predictive checks
  pp_check(b1, ndraws=100) + xlab('p(match)') + ggtitle('Posterior predictive density overlay');
  grid.arrange(pp_check(b1, type='stat', sta ='min') + xlab('p(match)') + ggtitle('Posterior predictive distribution of minimum values'),
               pp_check(b1, type='stat', stat='mean') + xlab('p(match)') + ggtitle('Posterior predictive distribution of means'),
               pp_check(b1, type='stat', stat='max') + xlab('p(match)') + ggtitle('Posterior predictive distribution of maximum values'),
               ncol=1); # seems fine (the observed, y, does fall in the predicted distributions, y_{rep} for the max, mean and min)
  
  # save the model:
  b1_res_web <- list("b_priors"=b_priors, "b1"=b1);
  saveRDS(b1_res_web, "./cached_results/b1_res_web.RData", compress="xz");
} else
{
  b1_res_web <- readRDS("./cached_results/b1_res_web.RData");
}

Priors As per various guidelines (see, for example, https://github.com/stan-dev/stan/wiki/Prior-Choice-Recommendations) and experiments, we will use student_t(5, 0, 2.5) for the slopes β, and lkj(2) for the random effects correlation structure, keeping the default brms priors for everything else, as they seem both sane and to perform well in our prior predictive checks. So, the priors we use are:

get_prior(b1_res_web$b_priors);

grid.arrange(pp_check(b1_res_web$b_priors, type='stat', sta ='min') + xlab('p(match)') + ggtitle('Prior predictive distribution of minimum values'),
             pp_check(b1_res_web$b_priors, type='stat', stat='mean') + xlab('p(match)') + ggtitle('Prior predictive distribution of means'),
             pp_check(b1_res_web$b_priors, type='stat', stat='max') + xlab('p(match)') + ggtitle('Prior predictive distribution of maximum values'),
             ncol=1); # seems fine but our value is a bit extreme

## Using all posterior draws for ppc type 'stat' by default.
## Using all posterior draws for ppc type 'stat' by default.
## Using all posterior draws for ppc type 'stat' by default.

## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.

Figure 3. Prior predictive checks. Figure generated using R version 4.3.3 (2024-02-29)

and it can be seen that they are quite ok, even if the observed p(match) is a bit extreme, but still within what the prior distribution covers.

Convergence The model converges well (high ESS and \(\hat{R} = 1.0\)):

summary(b1_res_web$b1);

##  Family: bernoulli 
##   Links: mu = logit 
## Formula: Match ~ 1 + Order + (1 + Order | Language) + (1 + Order | Family) + (1 + Order | Autotyp_Area) 
##    Data: web (Number of observations: 903) 
##   Draws: 4 chains, each with iter = 10000; warmup = 4000; thin = 2;
##          total post-warmup draws = 12000
## 
## Multilevel Hyperparameters:
## ~Autotyp_Area (Number of levels: 6) 
##                       Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS
## sd(Intercept)             0.43      0.38     0.02     1.40 1.00     7044
## sd(Order1)                0.68      0.62     0.02     2.31 1.00     8744
## cor(Intercept,Order1)     0.05      0.45    -0.80     0.84 1.00    10857
##                       Tail_ESS
## sd(Intercept)             9633
## sd(Order1)               10223
## cor(Intercept,Order1)    10571
## 
## ~Family (Number of levels: 9) 
##                       Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS
## sd(Intercept)             0.36      0.30     0.01     1.12 1.00     8623
## sd(Order1)                0.69      0.55     0.03     2.06 1.00     8960
## cor(Intercept,Order1)     0.03      0.45    -0.80     0.83 1.00    10376
##                       Tail_ESS
## sd(Intercept)             9633
## sd(Order1)                9252
## cor(Intercept,Order1)    10974
## 
## ~Language (Number of levels: 25) 
##                       Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS
## sd(Intercept)             0.55      0.21     0.16     1.02 1.00     6419
## sd(Order1)                0.93      0.36     0.26     1.70 1.00     7836
## cor(Intercept,Order1)     0.42      0.33    -0.31     0.92 1.00     8766
##                       Tail_ESS
## sd(Intercept)             5539
## sd(Order1)                6938
## cor(Intercept,Order1)     9511
## 
## Regression Coefficients:
##           Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## Intercept     2.03      0.36     1.28     2.73 1.00     9372     9214
## Order1       -1.04      0.61    -2.34     0.12 1.00    10071    10084
## 
## Draws were sampled using sampling(NUTS). For each parameter, Bulk_ESS
## and Tail_ESS are effective sample size measures, and Rhat is the potential
## scale reduction factor on split chains (at convergence, Rhat = 1).

mcmc_plot(b1_res_web$b1, type="trace");

## No divergences to plot.

Figure 4. Trace of the fitting process. Figure generated using R version 4.3.3 (2024-02-29)

Posterior predictive checks seems to be perfectly fine:

pp_check(b1_res_web$b1, ndraws=500) + xlab('p(match)') + ggtitle('Posterior predictive density overlay');

Figure 5. Posterior predictive density overlay for n=500 draws. Figure generated using R version 4.3.3 (2024-02-29)

grid.arrange(pp_check(b1_res_web$b1, type='stat', sta ='min') + xlab('p(match)') + ggtitle('Posterior predictive distribution of minimum values'),
             pp_check(b1_res_web$b1, type='stat', stat='mean') + xlab('p(match)') + ggtitle('Posterior predictive distribution of means'),
             pp_check(b1_res_web$b1, type='stat', stat='max') + xlab('p(match)') + ggtitle('Posterior predictive distribution of maximum values'),
             ncol=1); 

## Using all posterior draws for ppc type 'stat' by default.
## Using all posterior draws for ppc type 'stat' by default.
## Using all posterior draws for ppc type 'stat' by default.

## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.

Figure 6. Posterior predictive distribution. Figure generated using R version 4.3.3 (2024-02-29)

Fixed effects:

mcmc_plot(b1_res_web$b1);

Figure 7. Fixed effects’ posterior distributions. Figure generated using R version 4.3.3 (2024-02-29)

plot_model(b1_res_web$b1, type="emm", terms=c("Order"));

Figure 8. Fixed effects. Figure generated using R version 4.3.3 (2024-02-29)

Intercept The intercept is 2.03 95%CI [1.28, 2.73] on the log odds ratio (LOR) scale, which is clearly > 0; formal hypothesis testing against 0 is p(α>0)=1.00*; this translates into a probability of a match p(match)=88.4% with 95%CI [78.3%, 93.9%] ≫ 50%.

Order The slope of Order is -1.04 95%CI [-2.34, 0.12] on the log odds ratio (LOR) scale, which includes 0; formal hypothesis testing against 0 is p(β<0)=0.96*. In terms of probabilities, presenting [l] first (Order == “l_first”) results in a decrease by 15.5% [0.666%, 52.4%] in the probability of a match from the 88.4% [78.3%, 93.9%] when [r] is presented first 73.0% [38.3%, 91.9%] when presenting [r] first.

Random effects:

#get_variables(b1_res_web$b1);
grid.arrange(
  # 1 | Language:
  b1_res_web$b1 %>%
    spread_draws(b_Intercept, r_Language[Language, Order]) %>%
    filter(Order == "Intercept") %>%
    mutate(condition_mean = r_Language) %>%
    ggplot(aes(y = Language, x = condition_mean, fill = after_stat(x < 0))) +
    stat_halfeye(alpha=0.75, point_interval=median_qi, .width=c(.50, .95)) +
    scale_fill_manual(values = c("lightsalmon", "skyblue"), name=NULL, labels=c("≥ 0", "< 0")) +
    geom_vline(xintercept = 0, linetype = "dashed") +
    ggtitle("1 | Language") + xlab("Intercept") + ylab(NULL) + xlim(-2,2) +
    theme(legend.position="bottom") +
    NULL,
  
  # 1 | Family:
  b1_res_web$b1 %>%
    spread_draws(b_Intercept, r_Family[Family, Order]) %>%
    filter(Order == "Intercept") %>%
    mutate(condition_mean = r_Family) %>%
    ggplot(aes(y = Family, x = condition_mean, fill = after_stat(x < 0))) +
    stat_halfeye(alpha=0.75, point_interval=median_qi, .width=c(.50, .95)) +
    scale_fill_manual(values = c("lightsalmon", "skyblue"), name=NULL, labels=c("≥ 0", "< 0")) +
    geom_vline(xintercept = 0, linetype = "dashed") +
    ggtitle("1 | Family") + xlab("Intercept") + ylab(NULL) + xlim(-2,2) +
    theme(legend.position="bottom") +
    NULL,
  
  # 1 | Autotyp_Area:
  b1_res_web$b1 %>%
    spread_draws(b_Intercept, r_Autotyp_Area[Autotyp_Area, Order]) %>%
    filter(Order == "Intercept") %>%
    mutate(condition_mean = r_Autotyp_Area) %>%
    ggplot(aes(y = Autotyp_Area, x = condition_mean, fill = after_stat(x < 0))) +
    stat_halfeye(alpha=0.75, point_interval=median_qi, .width=c(.50, .95)) +
    scale_fill_manual(values = c("lightsalmon", "skyblue"), name=NULL, labels=c("≥ 0", "< 0")) +
    geom_vline(xintercept = 0, linetype = "dashed") +
    ggtitle("1 | Area") + xlab("Intercept") + ylab(NULL) + xlim(-2,2) +
    theme(legend.position="bottom") +
    NULL,
  
  # Order | Language:
  b1_res_web$b1 %>%
    spread_draws(b_Intercept, r_Language[Language, Order]) %>%
    filter(Order == "Order1") %>%
    mutate(condition_mean = r_Language) %>%
    ggplot(aes(y = Language, x = condition_mean, fill = after_stat(x < 0))) +
    stat_halfeye(alpha=0.75, point_interval=median_qi, .width=c(.50, .95)) +
    scale_fill_manual(values = c("lightsalmon", "skyblue"), name=NULL, labels=c("≥ 0", "< 0")) +
    geom_vline(xintercept = 0, linetype = "dashed") +
    ggtitle("Order | Language") + xlab("Slope") + ylab(NULL) + xlim(-4, 4) +
    theme(legend.position="bottom") +
    NULL,
  
  # Order | Family:
  b1_res_web$b1 %>%
    spread_draws(b_Intercept, r_Family[Family, Order]) %>%
    filter(Order == "Order1") %>%
    mutate(condition_mean = r_Family) %>%
    ggplot(aes(y = Family, x = condition_mean, fill = after_stat(x < 0))) +
    stat_halfeye(alpha=0.75, point_interval=median_qi, .width=c(.50, .95)) +
    scale_fill_manual(values = c("lightsalmon", "skyblue"), name=NULL, labels=c("≥ 0", "< 0")) +
    geom_vline(xintercept = 0, linetype = "dashed") +
    ggtitle("Order | Family") + xlab("Slope") + ylab(NULL) + xlim(-3, 3) + 
    theme(legend.position="bottom") +
    NULL,
  
  # Order | Autotyp_Area:
  b1_res_web$b1 %>%
    spread_draws(b_Intercept, r_Autotyp_Area[Autotyp_Area, Order]) %>%
    filter(Order == "Order1") %>%
    mutate(condition_mean = r_Autotyp_Area) %>%
    ggplot(aes(y = Autotyp_Area, x = condition_mean, fill = after_stat(x < 0))) +
    stat_halfeye(alpha=0.75, point_interval=median_qi, .width=c(.50, .95)) +
    scale_fill_manual(values = c("lightsalmon", "skyblue"), name=NULL, labels=c("≥ 0", "< 0")) +
    geom_vline(xintercept = 0, linetype = "dashed") +
    ggtitle("Order | Area") + xlab("Slope") + ylab(NULL) + xlim(-3, 3) +  
    theme(legend.position="bottom") +
    NULL,
  
  nrow=2);

Figure 9. Posterior estimaes of the random effects showing the median (dot), 50% (thick line) and 95% (thin lines) quantiles. Also showing 0 as the vertical dotted line and negative (blue) vs positive (red) values. Figure generated using R version 4.3.3 (2024-02-29)

Interpretation: The overall probability of a “perfect” match in the online experiment, while controlling for Language, Family and Area, is clearly much higher than the conservative change level of 50%, being estimated as about 88% with a 95%CI of about [78%, 94%], showing that there is a very strong tendency across our participants to associate [r] with the wavy line and [l] with the straight line. The order of presentation does seem to matter, in the sense that presenting [r] first increases the chance of a match by about 15% with a wide 95% [0.7%, 52%].

3.3.6 Model 2: do the characteristics of the languages matter?

if( !file.exists("./cached_results/b2_res_web.RData") )
{
  # the model:
  b2 <- brm(Match ~ 1 + Order + # random slope
              r_l_distinction_L1 + r_l_distinction_L2 + trill_real_L1 + # only random intercepts for these
              trill_real_L2 + # but random slope for this one
              (1 + Order + trill_real_L2 | Language) +
              (1 + Order + trill_real_L2 | Family) +
              (1 + Order + trill_real_L2 | Autotyp_Area),
            data=web,
            family=bernoulli(link='logit'),
            prior=c(brms::set_prior("student_t(5, 0, 2.5)", class="b"), 
                    brms::set_prior("lkj(2)", class="cor")),
            save_pars=save_pars(all=TRUE), # needed for Bayes factors
            sample_prior=TRUE,  # needed for hypotheses tests
            seed=998, cores=brms_ncores, iter=10000, warmup=4000, thin=2, control=list(adapt_delta=0.999, max_treedepth=13));
  summary(b2); mcmc_plot(b2, type="trace"); mcmc_plot(b2); # very decent
  mcmc_plot(b2, type="trace");
  bayestestR::hdi(b2, ci=0.95);
  hypothesis(b2, c("Order1 = 0", "r_l_distinction_L1distinct = 0", "r_l_distinction_L2distinct = 0", "trill_real_L1yes = 0", "trill_real_L2yes = 0")); 
  # posterior predictive checks
  pp_check(b2, ndraws=100) + xlab('p(match)') + ggtitle('Posterior predictive density overlay');
  grid.arrange(pp_check(b2, type='stat', sta ='min') + xlab('p(match)') + ggtitle('Posterior predictive distribution of minimum values') + xlim(0,1),
               pp_check(b2, type='stat', stat='mean') + xlab('p(match)') + ggtitle('Posterior predictive distribution of means') + xlim(0,1),
               pp_check(b2, type='stat', stat='max') + xlab('p(match)') + ggtitle('Posterior predictive distribution of maximum values') + xlim(0,1),
               ncol=1); # seems perfect
  
  # save the model:
  b2_res_web <- list("b2"=b2);
  saveRDS(b2_res_web, "./cached_results/b2_res_web.RData", compress="xz");
} else
{
  b2_res_web <- readRDS("./cached_results/b2_res_web.RData");
}

Priors are as before.

Convergence The model converges well (high ESS and \(\hat{R} = 1.0\)):

summary(b2_res_web$b2);

##  Family: bernoulli 
##   Links: mu = logit 
## Formula: Match ~ 1 + Order + r_l_distinction_L1 + r_l_distinction_L2 + trill_real_L1 + trill_real_L2 + (1 + Order + trill_real_L2 | Language) + (1 + Order + trill_real_L2 | Family) + (1 + Order + trill_real_L2 | Autotyp_Area) 
##    Data: web (Number of observations: 781) 
##   Draws: 4 chains, each with iter = 10000; warmup = 4000; thin = 2;
##          total post-warmup draws = 12000
## 
## Multilevel Hyperparameters:
## ~Autotyp_Area (Number of levels: 6) 
##                                 Estimate Est.Error l-95% CI u-95% CI Rhat
## sd(Intercept)                       0.60      0.55     0.02     2.00 1.00
## sd(Order1)                          1.02      0.88     0.04     3.28 1.00
## sd(trill_real_L2yes)                1.16      1.00     0.05     3.74 1.00
## cor(Intercept,Order1)               0.08      0.42    -0.72     0.81 1.00
## cor(Intercept,trill_real_L2yes)    -0.02      0.41    -0.77     0.74 1.00
## cor(Order1,trill_real_L2yes)       -0.00      0.42    -0.76     0.76 1.00
##                                 Bulk_ESS Tail_ESS
## sd(Intercept)                       7159     9073
## sd(Order1)                          6980     8396
## sd(trill_real_L2yes)                7534     8426
## cor(Intercept,Order1)              10077    10771
## cor(Intercept,trill_real_L2yes)    10390    10936
## cor(Order1,trill_real_L2yes)       10479    10735
## 
## ~Family (Number of levels: 9) 
##                                 Estimate Est.Error l-95% CI u-95% CI Rhat
## sd(Intercept)                       0.86      0.58     0.04     2.23 1.00
## sd(Order1)                          1.05      0.79     0.04     2.98 1.00
## sd(trill_real_L2yes)                0.91      0.85     0.03     3.13 1.00
## cor(Intercept,Order1)               0.14      0.41    -0.68     0.83 1.00
## cor(Intercept,trill_real_L2yes)    -0.08      0.41    -0.81     0.71 1.00
## cor(Order1,trill_real_L2yes)       -0.03      0.41    -0.78     0.75 1.00
##                                 Bulk_ESS Tail_ESS
## sd(Intercept)                       5893     7042
## sd(Order1)                          6711     7431
## sd(trill_real_L2yes)                9535     9792
## cor(Intercept,Order1)               9073     9221
## cor(Intercept,trill_real_L2yes)    10272     9746
## cor(Order1,trill_real_L2yes)       10743    10102
## 
## ~Language (Number of levels: 25) 
##                                 Estimate Est.Error l-95% CI u-95% CI Rhat
## sd(Intercept)                       0.39      0.27     0.02     1.01 1.00
## sd(Order1)                          0.91      0.46     0.09     1.88 1.00
## sd(trill_real_L2yes)                0.79      0.47     0.05     1.85 1.00
## cor(Intercept,Order1)               0.21      0.39    -0.63     0.85 1.00
## cor(Intercept,trill_real_L2yes)    -0.07      0.40    -0.77     0.71 1.00
## cor(Order1,trill_real_L2yes)        0.19      0.37    -0.58     0.82 1.00
##                                 Bulk_ESS Tail_ESS
## sd(Intercept)                       5511     7578
## sd(Order1)                          5891     4854
## sd(trill_real_L2yes)                5811     5807
## cor(Intercept,Order1)               5995     8490
## cor(Intercept,trill_real_L2yes)     8223     9753
## cor(Order1,trill_real_L2yes)        8951    10113
## 
## Regression Coefficients:
##                            Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS
## Intercept                      3.58      2.63    -0.67     9.41 1.00     9134
## Order1                        -0.87      0.80    -2.54     0.68 1.00     9786
## r_l_distinction_L1distinct     0.25      0.95    -1.71     2.05 1.00     9300
## r_l_distinction_L2distinct    -1.16      2.53    -6.80     2.83 1.00     8942
## trill_real_L1yes              -1.05      0.42    -1.91    -0.25 1.00    10194
## trill_real_L2yes               0.13      0.95    -1.61     2.26 1.00     8981
##                            Tail_ESS
## Intercept                      7826
## Order1                         9854
## r_l_distinction_L1distinct    10079
## r_l_distinction_L2distinct     7859
## trill_real_L1yes              10650
## trill_real_L2yes               8416
## 
## Draws were sampled using sampling(NUTS). For each parameter, Bulk_ESS
## and Tail_ESS are effective sample size measures, and Rhat is the potential
## scale reduction factor on split chains (at convergence, Rhat = 1).

mcmc_plot(b2_res_web$b2, type="trace");

## No divergences to plot.

Figure 10. Trace of the fitting process. Figure generated using R version 4.3.3 (2024-02-29)

Posterior predictive checks seems to be perfectly fine:

pp_check(b2_res_web$b2, ndraws=500) + xlab('p(match)') + ggtitle('Posterior predictive density overlay');

Figure 11. Posterior predictive density overlay for n=500 draws. Figure generated using R version 4.3.3 (2024-02-29)

grid.arrange(pp_check(b2_res_web$b2, type='stat', sta ='min') + xlab('p(match)') + ggtitle('Posterior predictive distribution of minimum values') + xlim(0,1),
             pp_check(b2_res_web$b2, type='stat', stat='mean') + xlab('p(match)') + ggtitle('Posterior predictive distribution of means') + xlim(0,1),
             pp_check(b2_res_web$b2, type='stat', stat='max') + xlab('p(match)') + ggtitle('Posterior predictive distribution of maximum values') + xlim(0,1),
             ncol=1); 

## Using all posterior draws for ppc type 'stat' by default.
## Using all posterior draws for ppc type 'stat' by default.
## Using all posterior draws for ppc type 'stat' by default.

## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.

Figure 12. Posterior predictive distribution (the x-axis scale was forced to [0,1]). Figure generated using R version 4.3.3 (2024-02-29)

Fixed effects:

mcmc_plot(b2_res_web$b2);

Figure 13. Fixed effects’ posterior distributions. Figure generated using R version 4.3.3 (2024-02-29)

grid.arrange(plot_model(b2_res_web$b2, type="emm", terms=c("Order")),
             plot_model(b2_res_web$b2, type="emm", terms=c("r_l_distinction_L1")),
             plot_model(b2_res_web$b2, type="emm", terms=c("r_l_distinction_L2")),
             plot_model(b2_res_web$b2, type="emm", terms=c("trill_real_L1")),
             plot_model(b2_res_web$b2, type="emm", terms=c("trill_real_L2")),
             ncol=2);

Figure 14. Fixed effects. Figure generated using R version 4.3.3 (2024-02-29)

Intercept The intercept is 3.58 95%CI [-0.67, 9.41] on the log odds ratio (LOR) scale, which is > 0 but has very large uncertainties so that the 95%CI includes 0; formal hypothesis testing against 0 is p(α>0)=0.95; this translates into a probability of a match p(match)=97.3% 95%CI [33.9%, 100.0%] which includes 50% in the 95%CI.

Order The slope of Order is -0.87 95%CI [-2.54, 0.68] on the log odds ratio (LOR) scale, which includes 0; formal hypothesis testing against 0 is p(β<0)=0.89. In terms of probabilities, presenting [l] first (Order1 == “l_first”) results in a decrease by 3.5% [0.004%, 30.0%] in the probability of a match from the 97.3% [33.9%, 100.0%] when [r] is presented first 93.8% [14.8%, 100.0%] when [l] is presented first.

Are [r] and [l] distinct in the L1(s) spoken by the participant? The slope of r_l_distinction_L1 is 0.25 95%CI [-1.71, 2.05] on the log odds ratio (LOR) scale, which includes 0; formal hypothesis testing against 0 is p(β>0)=0.62.

Are [r] and [l] distinct in the L2(s) spoken by the participant? The slope of r_l_distinction_L2 is -1.16 95%CI [-6.80, 2.83] on the log odds ratio (LOR) scale, which includes 0; formal hypothesis testing against 0 is p(β<0)=0.66.

Is [r] the main allophone in the L1(s) spoken by the participant? The slope of trill_real_L1 is -1.05 95%CI [-1.91, -0.25] on the log odds ratio (LOR) scale, which is clearly < 0; formal hypothesis testing against 0 is p(β<0)=0.99*. In terms of probabilities, having [r] as the main allophone (trill_real_L1 == “yes”) results in a decrease by 4.6% [0.002%, 26.8%] in the probability of a match from the 97.3% [33.9%, 100.0%] when [r] is not the primary allophone to 92.6% [14.6%, 100.0%] when it is.

Is [r] the main allophone in the L2(s) spoken by the participant? The slope of trill_real_L2 is 0.13 95%CI [-1.61, 2.26] on the log odds ratio (LOR) scale, which includes 0; formal hypothesis testing against 0 is p(β>0)=0.53.

Random effects:

#get_variables(b2_res_web$b2);
grid.arrange(
  # 1 | Language:
  b2_res_web$b2 %>%
    spread_draws(b_Intercept, r_Language[Language, Order]) %>%
    filter(Order == "Intercept") %>%
    mutate(condition_mean = r_Language) %>%
    ggplot(aes(y = Language, x = condition_mean, fill = after_stat(x < 0))) +
    stat_halfeye(alpha=0.75, point_interval=median_qi, .width=c(.50, .95)) +
    scale_fill_manual(values = c("lightsalmon", "skyblue"), name=NULL, labels=c("≥ 0", "< 0")) +
    geom_vline(xintercept = 0, linetype = "dashed") +
    ggtitle("1 | Language") + xlab("Intercept") + ylab(NULL) + xlim(-2,2) +
    theme(legend.position="bottom") +
    NULL,
  
  # 1 | Family:
  b2_res_web$b2 %>%
    spread_draws(b_Intercept, r_Family[Family, Order]) %>%
    filter(Order == "Intercept") %>%
    mutate(condition_mean = r_Family) %>%
    ggplot(aes(y = Family, x = condition_mean, fill = after_stat(x < 0))) +
    stat_halfeye(alpha=0.75, point_interval=median_qi, .width=c(.50, .95)) +
    scale_fill_manual(values = c("lightsalmon", "skyblue"), name=NULL, labels=c("≥ 0", "< 0")) +
    geom_vline(xintercept = 0, linetype = "dashed") +
    ggtitle("1 | Family") + xlab("Intercept") + ylab(NULL) + xlim(-2,2) +
    theme(legend.position="bottom") +
    NULL,
  
  # 1 | Autotyp_Area:
  b2_res_web$b2 %>%
    spread_draws(b_Intercept, r_Autotyp_Area[Autotyp_Area, Order]) %>%
    filter(Order == "Intercept") %>%
    mutate(condition_mean = r_Autotyp_Area) %>%
    ggplot(aes(y = Autotyp_Area, x = condition_mean, fill = after_stat(x < 0))) +
    stat_halfeye(alpha=0.75, point_interval=median_qi, .width=c(.50, .95)) +
    scale_fill_manual(values = c("lightsalmon", "skyblue"), name=NULL, labels=c("≥ 0", "< 0")) +
    geom_vline(xintercept = 0, linetype = "dashed") +
    ggtitle("1 | Area") + xlab("Intercept") + ylab(NULL) + xlim(-2,2) +
    theme(legend.position="bottom") +
    NULL,
  
  # Order | Language:
  b2_res_web$b2 %>%
    spread_draws(b_Intercept, r_Language[Language, Order]) %>%
    filter(Order == "Order1") %>%
    mutate(condition_mean = r_Language) %>%
    ggplot(aes(y = Language, x = condition_mean, fill = after_stat(x < 0))) +
    stat_halfeye(alpha=0.75, point_interval=median_qi, .width=c(.50, .95)) +
    scale_fill_manual(values = c("lightsalmon", "skyblue"), name=NULL, labels=c("≥ 0", "< 0")) +
    geom_vline(xintercept = 0, linetype = "dashed") +
    ggtitle("Order | Language") + xlab("Slope") + ylab(NULL) + xlim(-4, 4) +
    theme(legend.position="bottom") +
    NULL,
  
  # Order | Family:
  b2_res_web$b2 %>%
    spread_draws(b_Intercept, r_Family[Family, Order]) %>%
    filter(Order == "Order1") %>%
    mutate(condition_mean = r_Family) %>%
    ggplot(aes(y = Family, x = condition_mean, fill = after_stat(x < 0))) +
    stat_halfeye(alpha=0.75, point_interval=median_qi, .width=c(.50, .95)) +
    scale_fill_manual(values = c("lightsalmon", "skyblue"), name=NULL, labels=c("≥ 0", "< 0")) +
    geom_vline(xintercept = 0, linetype = "dashed") +
    ggtitle("Order | Family") + xlab("Slope") + ylab(NULL) + xlim(-3, 3) + 
    theme(legend.position="bottom") +
    NULL,
  
  # Order | Autotyp_Area:
  b2_res_web$b2 %>%
    spread_draws(b_Intercept, r_Autotyp_Area[Autotyp_Area, Order]) %>%
    filter(Order == "Order1") %>%
    mutate(condition_mean = r_Autotyp_Area) %>%
    ggplot(aes(y = Autotyp_Area, x = condition_mean, fill = after_stat(x < 0))) +
    stat_halfeye(alpha=0.75, point_interval=median_qi, .width=c(.50, .95)) +
    scale_fill_manual(values = c("lightsalmon", "skyblue"), name=NULL, labels=c("≥ 0", "< 0")) +
    geom_vline(xintercept = 0, linetype = "dashed") +
    ggtitle("Order | Area") + xlab("Slope") + ylab(NULL) + xlim(-3, 3) +  
    theme(legend.position="bottom") +
    NULL,
  
  # trill_real_L2 | Language:
  b2_res_web$b2 %>%
    spread_draws(b_Intercept, r_Language[Language, trill_real_L2]) %>%
    filter(trill_real_L2 == "trill_real_L2yes") %>%
    mutate(condition_mean = r_Language) %>%
    ggplot(aes(y = Language, x = condition_mean, fill = after_stat(x < 0))) +
    stat_halfeye(alpha=0.75, point_interval=median_qi, .width=c(.50, .95)) +
    scale_fill_manual(values = c("lightsalmon", "skyblue"), name=NULL, labels=c("≥ 0", "< 0")) +
    geom_vline(xintercept = 0, linetype = "dashed") +
    ggtitle("trill_real_L2 | Language") + xlab("Slope") + ylab(NULL) + xlim(-4, 4) +
    theme(legend.position="bottom") +
    NULL,
  
  # trill_real_L2 | Family:
  b2_res_web$b2 %>%
    spread_draws(b_Intercept, r_Family[Family, trill_real_L2]) %>%
    filter(trill_real_L2 == "trill_real_L2yes") %>%
    mutate(condition_mean = r_Family) %>%
    ggplot(aes(y = Family, x = condition_mean, fill = after_stat(x < 0))) +
    stat_halfeye(alpha=0.75, point_interval=median_qi, .width=c(.50, .95)) +
    scale_fill_manual(values = c("lightsalmon", "skyblue"), name=NULL, labels=c("≥ 0", "< 0")) +
    geom_vline(xintercept = 0, linetype = "dashed") +
    ggtitle("trill_real_L2 | Family") + xlab("Slope") + ylab(NULL) + xlim(-3, 3) + 
    theme(legend.position="bottom") +
    NULL,
  
  # trill_real_L2 | Autotyp_Area:
  b2_res_web$b2 %>%
    spread_draws(b_Intercept, r_Autotyp_Area[Autotyp_Area, trill_real_L2]) %>%
    filter(trill_real_L2 == "trill_real_L2yes") %>%
    mutate(condition_mean = r_Autotyp_Area) %>%
    ggplot(aes(y = Autotyp_Area, x = condition_mean, fill = after_stat(x < 0))) +
    stat_halfeye(alpha=0.75, point_interval=median_qi, .width=c(.50, .95)) +
    scale_fill_manual(values = c("lightsalmon", "skyblue"), name=NULL, labels=c("≥ 0", "< 0")) +
    geom_vline(xintercept = 0, linetype = "dashed") +
    ggtitle("trill_real_L2 | Area") + xlab("Slope") + ylab(NULL) + xlim(-3, 3) +  
    theme(legend.position="bottom") +
    NULL,
  
  nrow=3);

Figure 15. Posterior estimaes of the random effects showing the median (dot), 50% (thick line) and 95% (thin lines) quantiles. Also showing 0 as the vertical dotted line and negative (blue) vs positive (red) values. Figure generated using R version 4.3.3 (2024-02-29)

Interpretation: The overall probability of a “perfect” match in the online experiment when adding all the relevant [r]-related predictors for the L1(s) and L2(s) known to the participants, while controlling for Language, Family and Area, seem to still be higher than the conservative change level of 50%, being estimated as about 97% but with a much wider 95%CI of about [34%, 100.0%] which now does include 50%. The only predictor that seems to make a difference is if [r] is the main allophone in the L1(s) known by the participant, in that those participants that have an L1 with [r] as the main allophone have a slightly smaller probability of a perfect match by about 4.5% [0%, 27%] to about 93% [15%, 100%] than the others. Here, order does not seem to make clear a difference, but there is still a suggestion that presenting [r] first increases the probability of a match by about 3.5% [0%, 30%].

4 Field experiment

4.1 Descriptive statistics

First, how many participants?

nrow(field)

## [1] 127

Sex division

table(field$Sex)

## 
##  f  m 
## 91 36

Ages

summary(field$Age)

##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##   18.00   19.00   20.00   28.55   34.50   75.00

First, how many languages?

field %>% count(Language) %>% nrow()

## [1] 6

Does this number correspond with the L1s?

field %>% count(Name) %>% nrow()

## [1] 6

How many families?

field %>% count(Family) %>% nrow()

## [1] 4

How many have the R/L distinction in the L1 among the languages?

field %>% count(Name, r_l_distinction_L1) %>% count(r_l_distinction_L1)

How many really use the alveolar trill in L1 among the languages?

field %>% count(Name, trill_real_L1) %>% count(trill_real_L1)

How many really have the alveolar trill in L1 as an allophone among the languages?

field %>% count(Name, trill_occ_L1) %>% count(trill_occ_L1)

What about the same questions for L2. But this will not neatly sum up to 25, due to various possible scenarios for L2 within a specific L1.

How many have the R/L distinction in the L2 among the languages?

field %>% count(Name, r_l_distinction_L2) %>% count(r_l_distinction_L2)

How many really use the alveolar trill in L2 among the languages?

field %>% count(Name, trill_real_L2) %>% count(trill_real_L2)

How many really have the alveolar trill in L2 as an allophone among the languages?

field %>% count(Name, trill_occ_L2) %>% count(trill_occ_L2)

What is the grand average congruent behavior?

mean(field$Match)

## [1] 0.976378

97%!!!

What about only among those who have L1 without the distinction?

field %>%
  filter(r_l_distinction_L1 == "0") %>%
  summarize(mean_match = mean(Match, na.rm = TRUE))

100%! WOW.

What about only among those who have L1 without the distinction and no L2 that distinguishes?

field %>%
  filter(r_l_distinction_L1 == "0") %>%
  filter(!EnglishL2YesNo) %>% 
  filter(r_l_distinction_L2 == '0') %>% 
  summarize(mean_match = mean(Match, na.rm = TRUE))

There are no such people.

Compute average matching behavior per language and sort:

field_avg <- field %>%
  group_by(Language) %>% 
  summarize(M = mean(Match)) %>% 
  arrange(desc(M)) %>% 
  mutate(percent = round(M, 2) * 100,
         percent = str_c(percent, '%'))

# Show:

field_avg %>% print(n = Inf)

## # A tibble: 6 × 3
##   Language     M percent
##   <chr>    <dbl> <chr>  
## 1 BR       1     100%   
## 2 PA       1     100%   
## 3 VA       1     100%   
## 4 BE       0.982 98%    
## 5 DE       0.947 95%    
## 6 SR       0.923 92%

Check some demographics, also to report in the paper. First, the number of participants per language:

field %>% 
  count(Name, sort = TRUE) %>% print(n = Inf)

## # A tibble: 6 × 2
##   Name                     n
##   <chr>                <int>
## 1 English UK              55
## 2 Tashlhiyt Berber        20
## 3 German                  19
## 4 Brazilian Portuguese    13
## 5 Daakie                  12
## 6 Palikur                  8

Then, the number of L1 speakers who have R/L distinction vs. who don’t:

field %>% count(r_l_distinction_L1) %>%
  mutate(prop = n / sum(n),
         percent = round(prop, 2) * 100,
         percent = str_c(percent, '%'))

How many people do not have any L2?

# raw count of no L2; raw count of all ppl
sum(is.na(field$L2)); nrow(field)

## [1] 52

## [1] 127

# percentage no L2
sum(is.na(field$L2)) / nrow(field)

## [1] 0.4094488

# percentage with L2
1 - sum(is.na(field$L2)) / nrow(field)

## [1] 0.5905512

Check how many people knew English as their L2:

field %>% count(EnglishL2YesNo) %>% 
  mutate(prop = n / sum(n),
         percent = round(prop, 2) * 100,
         percent = str_c(percent, '%'))

Of those that don’t use a R/L distinction in their L1, how many use R/L distinction in their L2? (double-check if logic alright!)

field %>%
  filter(r_l_distinction_L1 == '0') %>% 
  count(r_l_distinction_L2 == '1') %>% 
  mutate(prop = n / sum(n),
         percent = round(prop, 2) * 100,
         percent = str_c(percent, '%'))

How many “pure” speakers were there?, i.e., those people that 1) don’t know English, 2) don’t use an L1 with a R/L distinction, and 3) don’t know an L2 that distinguishes R/L.

field %>% 
  filter(r_l_distinction_L1 == '0') %>%  # excludes English as well
  filter(!EnglishL2YesNo) %>% 
  filter(r_l_distinction_L2 == '0') %>% 
  nrow()

## [1] 0

None.

Now, some relevant plots:

Distribution of matches in the field data by Language.

Distribution of matches in the field data by Family.

Distribution of matches in the field data by Area.

Distribution of matches in the field data by [r]/[l] distinction in L1.

Distribution of matches in the field data by [r] is main allophone in L1.

Distribution of matches in the field data by [r]/[l] distinction in L2.

Distribution of matches in the field data by [r] is main allophone in L2.

4.2 Regression models

4.2.1 What predictors to include?

We use the same ideas as for the web experiment above, except that now Order does not exist anymore.

# Make sure we code the factors with the intended baseline levels and contrasts:
field$r_l_distinction_L1 <- factor(c("same", "distinct")[field$r_l_distinction_L1 + 1], levels=c("same", "distinct"));
field$trill_real_L1      <- factor(c("no", "yes")[field$trill_real_L1 + 1], levels=c("no", "yes"));
field$trill_occ_L1       <- factor(c("no", "yes")[field$trill_occ_L1 + 1], levels=c("no", "yes"));
field$r_l_distinction_L2 <- factor(c("same", "distinct")[field$r_l_distinction_L2 + 1], levels=c("same", "distinct"));
field$trill_real_L2      <- factor(c("no", "yes")[field$trill_real_L2 + 1], levels=c("no", "yes"));
field$trill_occ_L2       <- factor(c("no", "yes")[field$trill_occ_L2 + 1], levels=c("no", "yes"));

field %>% count(Match) %>% mutate(prop = paste0(round(100* n / sum(n),1),"%")) 

# match is quite unbalanced: 2.4% vs 97.6%

field %>% count(r_l_distinction_L1) %>% mutate(prop = paste0(round(100* n / sum(n),1),"%")) 

# highly imbalanced: 6.3% vs 93.7%

field %>% count(trill_real_L1) %>% mutate(prop = paste0(round(100* n / sum(n),1),"%")) 

# imbalanced: 74.8% vs 25.2%

field %>% count(trill_occ_L1) %>% mutate(prop = paste0(round(100* n / sum(n),1),"%")) 

# imbalanced: 6.3% vs 93.7%

## And for L2, just in case
field %>% count(r_l_distinction_L2) %>% mutate(prop = paste0(round(100* n / sum(n),1),"%")) 

# 41.7% missing, the rest: 0.8% vs 57.5%

field %>% count(trill_real_L2) %>% mutate(prop = paste0(round(100* n / sum(n),1),"%")) 

# 41.7% missing, the rest: 22% vs 36.2%

field %>% count(trill_occ_L2) %>% mutate(prop = paste0(round(100* n / sum(n),1),"%")) 

# 41.7% missing, the rest are all (58.3%) are "yes" ---> exclude

It can be seen that, in our data, trill_occ_L2 is 100% “yes” (for the non-missing data cases), which means that we will not include it in any further analyses.

However, r_l_distinction_L1 and trill_occ_L1 are perfectly correlated:

table(field$r_l_distinction_L1, field$trill_occ_L1);

##           
##             no yes
##   same       8   0
##   distinct   0 119

so we will only consider r_l_distinction_L1.

4.2.2 What random structure to use?

As for the web experiment above:

unique(field[,c("Name", "Family", "Autotyp_Area")]) %>% arrange(Autotyp_Area, Family, Name)

It seems that including Family and Area do not make any sense in this case.

r/l distinction in L1 (r_l_distinction_L1):

table(field$r_l_distinction_L1, field$Language)

##           
##            BE BR DE PA SR VA
##   same      0  0  0  8  0  0
##   distinct 55 20 19  0 13 12

table(field$r_l_distinction_L1, field$Family)

##           
##            Afro-Asiatic Arawakan Austronesian IE
##   same                0        8            0  0
##   distinct           20        0           12 87

table(field$r_l_distinction_L1, field$Autotyp_Area)

##           
##            Europe N Africa NE South America Oceania
##   same          0        0                8       0
##   distinct     87       20                0      12

This is perfectly uniform within the levels of each of the three factors, so random slopes are arguably not justified.

r/l distinction in L2 (r_l_distinction_L2):

table(field$r_l_distinction_L2, field$Language)

##           
##            BE BR DE PA SR VA
##   same      1  0  0  0  0  0
##   distinct 18 20 14  8  1 12

table(field$r_l_distinction_L2, field$Family)

##           
##            Afro-Asiatic Arawakan Austronesian IE
##   same                0        0            0  1
##   distinct           20        8           12 33

table(field$r_l_distinction_L2, field$Autotyp_Area)

##           
##            Europe N Africa NE South America Oceania
##   same          1        0                0       0
##   distinct     33       20                8      12

There is almost no “same” at all, so random slopes are arguably not justified.

[r] is main allophone in L1 (trill_real_L1):

table(field$trill_real_L1, field$Language)

##      
##       BE BR DE PA SR VA
##   no  55  0 19  8 13  0
##   yes  0 20  0  0  0 12

table(field$trill_real_L1, field$Family)

##      
##       Afro-Asiatic Arawakan Austronesian IE
##   no             0        8            0 87
##   yes           20        0           12  0

table(field$trill_real_L1, field$Autotyp_Area)

##      
##       Europe N Africa NE South America Oceania
##   no      87        0                8       0
##   yes      0       20                0      12

This is perfectly uniform within the levels of each of the three factors, so random slopes are arguably not justified.

[r] is main allophone in L2 (trill_real_L2):

table(field$trill_real_L2, field$Language)

##      
##       BE BR DE PA SR VA
##   no   9  0 10  8  1  0
##   yes 10 20  4  0  0 12

table(field$trill_real_L2, field$Family)

##      
##       Afro-Asiatic Arawakan Austronesian IE
##   no             0        8            0 20
##   yes           20        0           12 14

table(field$trill_real_L2, field$Autotyp_Area)

##      
##       Europe N Africa NE South America Oceania
##   no      20        0                8       0
##   yes     14       20                0      12

There is variation only within Indo-European/Europe, so random slopes are arguably not justified.

[r] is present in L1 (trill_occ_L1):

table(field$trill_occ_L1, field$Language)

##      
##       BE BR DE PA SR VA
##   no   0  0  0  8  0  0
##   yes 55 20 19  0 13 12

table(field$trill_occ_L1, field$Family)

##      
##       Afro-Asiatic Arawakan Austronesian IE
##   no             0        8            0  0
##   yes           20        0           12 87

table(field$trill_occ_L1, field$Autotyp_Area)

##      
##       Europe N Africa NE South America Oceania
##   no       0        0                8       0
##   yes     87       20                0      12

This is perfectly uniform within the levels of each of the three factors, so random slopes are arguably not justified.

Given these, we only include Language as a random effect, and we do not model any random slopes here.

4.2.3 Model 1: what is the overall probability of a match?

As for the web experiment, but this time without Order (so, the intercept-only model):

if( !file.exists("./cached_results/b1_res_field.RData") )
{
  # prior predictive checks:
  # what priors we can set:
  get_prior(Match ~ 1 +
              (1 | Language),
            data=field,
            family=bernoulli(link='logit'));
  b_priors <- brms::brm(Match ~ 1 +
                          (1 | Language),
                        data=field,
                        family=bernoulli(link='logit'),
                        sample_prior='only',  # needed for prior predictive checks
                        seed=998, cores=brms_ncores, iter=10000, warmup=4000, thin=2, control=list(adapt_delta=0.999, max_treedepth=13));
  grid.arrange(pp_check(b_priors, type='stat', sta ='min') + xlab('p(match)') + ggtitle('Prior predictive distribution of minimum values'),
               pp_check(b_priors, type='stat', stat='mean') + xlab('p(match)') + ggtitle('Prior predictive distribution of means'),
               pp_check(b_priors, type='stat', stat='max') + xlab('p(match)') + ggtitle('Prior predictive distribution of maximum values'),
               ncol=1); # seems fine but our value is quite extreme
  
  # the model:
  b1 <- brm(Match ~ 1 +
              (1 | Language),
            data=field,
            family=bernoulli(link='logit'),
            save_pars=save_pars(all=TRUE), # needed for Bayes factors
            sample_prior=TRUE,  # needed for hypotheses tests
            seed=998, cores=brms_ncores, iter=10000, warmup=4000, thin=2, control=list(adapt_delta=0.999, max_treedepth=13));
  summary(b1); mcmc_plot(b1, type="trace"); mcmc_plot(b1); # very decent
  mcmc_plot(b1, type="trace");
  bayestestR::hdi(b1, ci=0.95);
  # posterior predictive checks
  pp_check(b1, ndraws=100) + xlab('p(match)') + ggtitle('Posterior predictive density overlay');
  grid.arrange(pp_check(b1, type='stat', sta ='min') + xlab('p(match)') + ggtitle('Posterior predictive distribution of minimum values'),
               pp_check(b1, type='stat', stat='mean') + xlab('p(match)') + ggtitle('Posterior predictive distribution of means'),
               pp_check(b1, type='stat', stat='max') + xlab('p(match)') + ggtitle('Posterior predictive distribution of maximum values'),
               ncol=1); # far from ideal but seems ok...
  
  # save the model:
  b1_res_field <- list("b_priors"=b_priors, "b1"=b1);
  saveRDS(b1_res_field, "./cached_results/b1_res_field.RData", compress="xz");
} else
{
  b1_res_field <- readRDS("./cached_results/b1_res_field.RData");
}

grid.arrange(pp_check(b1_res_field$b_priors, type='stat', sta ='min') + xlab('p(match)') + ggtitle('Prior predictive distribution of minimum values'),
             pp_check(b1_res_field$b_priors, type='stat', stat='mean') + xlab('p(match)') + ggtitle('Prior predictive distribution of means'),
             pp_check(b1_res_field$b_priors, type='stat', stat='max') + xlab('p(match)') + ggtitle('Prior predictive distribution of maximum values'),
             ncol=1); # seems fine but our value is a bit extreme

## Using all posterior draws for ppc type 'stat' by default.
## Using all posterior draws for ppc type 'stat' by default.
## Using all posterior draws for ppc type 'stat' by default.

## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.

Figure 16. Prior predictive checks. Figure generated using R version 4.3.3 (2024-02-29)

and it can be seen that they are quite ok, even if the observed p(match) is quite extreme, but still within what the prior distribution covers.

Convergence The model converges well (high ESS and \(\hat{R} = 1.0\)):

summary(b1_res_field$b1);

##  Family: bernoulli 
##   Links: mu = logit 
## Formula: Match ~ 1 + (1 | Language) 
##    Data: field (Number of observations: 127) 
##   Draws: 4 chains, each with iter = 10000; warmup = 4000; thin = 2;
##          total post-warmup draws = 12000
## 
## Multilevel Hyperparameters:
## ~Language (Number of levels: 6) 
##               Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## sd(Intercept)     0.98      0.92     0.03     3.42 1.00     7041     7465
## 
## Regression Coefficients:
##           Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## Intercept     3.84      0.84     2.37     5.72 1.00     7664     7214
## 
## Draws were sampled using sampling(NUTS). For each parameter, Bulk_ESS
## and Tail_ESS are effective sample size measures, and Rhat is the potential
## scale reduction factor on split chains (at convergence, Rhat = 1).

mcmc_plot(b1_res_field$b1, type="trace");

## No divergences to plot.

Figure 17. Trace of the fitting process. Figure generated using R version 4.3.3 (2024-02-29)

Posterior predictive checks seems to be acceptable:

pp_check(b1_res_field$b1, ndraws=500) + xlab('p(match)') + ggtitle('Posterior predictive density overlay');

Figure 18. Posterior predictive density overlay for n=500 draws. Figure generated using R version 4.3.3 (2024-02-29)

grid.arrange(pp_check(b1_res_field$b1, type='stat', sta ='min') + xlab('p(match)') + ggtitle('Posterior predictive distribution of minimum values') + xlim(0,1),
             pp_check(b1_res_field$b1, type='stat', stat='mean') + xlab('p(match)') + ggtitle('Posterior predictive distribution of means') + xlim(0,1),
             pp_check(b1_res_field$b1, type='stat', stat='max') + xlab('p(match)') + ggtitle('Posterior predictive distribution of maximum values') + xlim(0,1),
             ncol=1); 

## Using all posterior draws for ppc type 'stat' by default.
## Using all posterior draws for ppc type 'stat' by default.
## Using all posterior draws for ppc type 'stat' by default.

## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.

Figure 19. Posterior predictive distribution. Figure generated using R version 4.3.3 (2024-02-29)

Fixed effects:

mcmc_plot(b1_res_field$b1);

Figure 20. Fixed effects’ posterior distributions. Figure generated using R version 4.3.3 (2024-02-29)

Intercept The intercept is 3.84 95%CI [2.37, 5.72] on the log odds ratio (LOR) scale, which is clearly » 0; formal hypothesis testing against 0 is p(α>0)=1.00*; this translates into a probability of a match p(match)=97.9% 95%CI [91.4%, 99.7%] » 50%.

Random effects:

#get_variables(b1_res_field$b1);
grid.arrange(
  # 1 | Language:
  b1_res_field$b1 %>%
    spread_draws(b_Intercept, r_Language[Language, Order]) %>%
    filter(Order == "Intercept") %>%
    mutate(condition_mean = r_Language) %>%
    ggplot(aes(y = Language, x = condition_mean, fill = after_stat(x < 0))) +
    stat_halfeye(alpha=0.75, point_interval=median_qi, .width=c(.50, .95)) +
    scale_fill_manual(values = c("lightsalmon", "skyblue"), name=NULL, labels=c("≥ 0", "< 0")) +
    geom_vline(xintercept = 0, linetype = "dashed") +
    ggtitle("1 | Language") + xlab("Intercept") + ylab(NULL) + xlim(-2,2) +
    theme(legend.position="bottom") +
    NULL,

  nrow=1);

Figure 21. Posterior estimaes of the random effects showing the median (dot), 50% (thick line) and 95% (thin lines) quantiles. Also showing 0 as the vertical dotted line and negative (blue) vs positive (red) values. Figure generated using R version 4.3.3 (2024-02-29)

Interpretation: The overall probability of a “perfect” match in the field experiment, while controlling for Language, is much higher than the conservative change level of 50%, being estimated as about 98%, with a 95% credible interval of about [91%, 100%] that excludes 50%, suggesting that there is a very strong tendency across our participants to associate [r] with the wavy line and [l] with the straight line.

4.2.4 Model 2: do the characteristics of the languages matter?

if( !file.exists("./cached_results/b2_res_field.RData") )
{
  # the model:
  b2 <- brm(Match ~ 1 + 
              r_l_distinction_L1 + r_l_distinction_L2 + 
              trill_real_L1 + trill_real_L2 + 
              (1 | Language),
            data=field,
            family=bernoulli(link='logit'),
            prior=c(brms::set_prior("student_t(5, 0, 2.5)", class="b")), # same slope prior as for web
            save_pars=save_pars(all=TRUE), # needed for Bayes factors
            sample_prior=TRUE,  # needed for hypotheses tests
            seed=998, cores=brms_ncores, iter=10000, warmup=4000, thin=2, control=list(adapt_delta=0.999, max_treedepth=13));
  summary(b2); mcmc_plot(b2, type="trace"); mcmc_plot(b2); # very decent
  mcmc_plot(b2, type="trace");
  bayestestR::hdi(b2, ci=0.95);
  hypothesis(b2, c("r_l_distinction_L1distinct = 0", "r_l_distinction_L2distinct = 0", "trill_real_L1yes = 0", "trill_real_L2yes = 0")); 
  # posterior predictive checks
  pp_check(b2, ndraws=100) + xlab('p(match)') + ggtitle('Posterior predictive density overlay');
  pp_check(b2, ndraws=100, prefix="ppd") + xlab('p(match)') + ggtitle('Posterior predictive density overlay without the observed data');
  grid.arrange(pp_check(b2, type='stat', sta ='min') + xlab('p(match)') + ggtitle('Posterior predictive distribution of minimum values') + xlim(0,1),
               pp_check(b2, type='stat', stat='mean') + xlab('p(match)') + ggtitle('Posterior predictive distribution of means') + xlim(0,1),
               pp_check(b2, type='stat', stat='max') + xlab('p(match)') + ggtitle('Posterior predictive distribution of maximum values') + xlim(0,1),
               ncol=1); # seems perfect
  
  # save the model:
  b2_res_field <- list("b2"=b2);
  saveRDS(b2_res_field, "./cached_results/b2_res_field.RData", compress="xz");
} else
{
  b2_res_field <- readRDS("./cached_results/b2_res_field.RData");
}

Priors are as before.

Convergence The model converges well (high ESS and \(\hat{R} = 1.0\)):

summary(b2_res_field$b2);

##  Family: bernoulli 
##   Links: mu = logit 
## Formula: Match ~ 1 + r_l_distinction_L1 + r_l_distinction_L2 + trill_real_L1 + trill_real_L2 + (1 | Language) 
##    Data: field (Number of observations: 74) 
##   Draws: 4 chains, each with iter = 10000; warmup = 4000; thin = 2;
##          total post-warmup draws = 12000
## 
## Multilevel Hyperparameters:
## ~Language (Number of levels: 6) 
##               Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## sd(Intercept)     2.05      2.38     0.06     7.30 1.00     9933     8547
## 
## Regression Coefficients:
##                            Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS
## Intercept                     10.26      7.57     0.17    27.28 1.00     7762
## r_l_distinction_L1distinct    -0.21      2.87    -6.14     5.30 1.00    10425
## r_l_distinction_L2distinct    -0.13      3.01    -6.39     5.68 1.00    11392
## trill_real_L1yes              -0.13      2.67    -5.49     5.16 1.00    11199
## trill_real_L2yes              -0.18      2.65    -5.55     4.98 1.00    10980
##                            Tail_ESS
## Intercept                      5939
## r_l_distinction_L1distinct     9293
## r_l_distinction_L2distinct     9593
## trill_real_L1yes               9418
## trill_real_L2yes               9461
## 
## Draws were sampled using sampling(NUTS). For each parameter, Bulk_ESS
## and Tail_ESS are effective sample size measures, and Rhat is the potential
## scale reduction factor on split chains (at convergence, Rhat = 1).

mcmc_plot(b2_res_field$b2, type="trace");

## No divergences to plot.

Figure 22. Trace of the fitting process. Figure generated using R version 4.3.3 (2024-02-29)

Posterior predictive checks are actually pretty good (please note that there are some curves hidden by the observed data, visible in the second plot):

pp_check(b2_res_field$b2, ndraws=500) + xlab('p(match)') + ggtitle('Posterior predictive density overlay');

Figure 23. Posterior predictive density overlay for n=500 draws. Figure generated using R version 4.3.3 (2024-02-29)

pp_check(b2_res_field$b2, ndraws=500, prefix="ppd") + xlab('p(match)') + ggtitle('Posterior predictive density overlay');

Figure 24. Posterior predictive density overlay without showinf the observed data for n=500 draws. Figure generated using R version 4.3.3 (2024-02-29)

grid.arrange(pp_check(b2_res_field$b2, type='stat', sta ='min') + xlab('p(match)') + ggtitle('Posterior predictive distribution of minimum values') + xlim(0,1),
             pp_check(b2_res_field$b2, type='stat', stat='mean') + xlab('p(match)') + ggtitle('Posterior predictive distribution of means') + xlim(0,1),
             pp_check(b2_res_field$b2, type='stat', stat='max') + xlab('p(match)') + ggtitle('Posterior predictive distribution of maximum values') + xlim(0,1),
             ncol=1); 

## Using all posterior draws for ppc type 'stat' by default.
## Using all posterior draws for ppc type 'stat' by default.
## Using all posterior draws for ppc type 'stat' by default.

## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.

Figure 25. Posterior predictive distribution (the x-axis scale was forced to [0,1]). Figure generated using R version 4.3.3 (2024-02-29)

Fixed effects:

mcmc_plot(b2_res_field$b2);

Figure 26. Fixed effects’ posterior distributions. Figure generated using R version 4.3.3 (2024-02-29)

grid.arrange(plot_model(b2_res_field$b2, type="emm", terms=c("r_l_distinction_L1")),
             plot_model(b2_res_field$b2, type="emm", terms=c("r_l_distinction_L2")),
             plot_model(b2_res_field$b2, type="emm", terms=c("trill_real_L1")),
             plot_model(b2_res_field$b2, type="emm", terms=c("trill_real_L2")),
             ncol=2);

Figure 27. Fixed effects. Figure generated using R version 4.3.3 (2024-02-29)

Intercept The intercept is 10.26 95%CI [0.17, 27.28] on the log odds ratio (LOR) scale, which is > 0 and with a 95%CI that does not includes 0; formal hypothesis testing against 0 is p(α>0)=0.98*; this translates into a probability of a match p(match)=100.0% 95%CI [54.1%, 100.0%] which does not include 50%.

Are [r] and [l] distinct in the L1(s) spoken by the participant? The slope of r_l_distinction_L1 is -0.21 95%CI [-6.14, 5.30] on the log odds ratio (LOR) scale, which includes 0; formal hypothesis testing against 0 is p(β<0)=0.52.

Are [r] and [l] distinct in the L2(s) spoken by the participant? The slope of r_l_distinction_L2 is -0.13 95%CI [-6.39, 5.68] on the log odds ratio (LOR) scale, which includes 0; formal hypothesis testing against 0 is p(β<0)=0.51.

Is [r] the main allophone in the L1(s) spoken by the participant? The slope of trill_real_L1 is -0.13 95%CI [-5.49, 5.16] on the log odds ratio (LOR) scale, which includes 0; formal hypothesis testing against 0 is p(β<0)=0.52.

Is [r] the main allophone in the L2(s) spoken by the participant? The slope of trill_real_L2 is -0.18 95%CI [-5.55, 4.98] on the log odds ratio (LOR) scale, which includes 0; formal hypothesis testing against 0 is p(β<0)=0.53.

Random effects:

#get_variables(b2_res_field$b2);
grid.arrange(
  # 1 | Language:
  b2_res_field$b2 %>%
    spread_draws(b_Intercept, r_Language[Language, Order]) %>%
    filter(Order == "Intercept") %>%
    mutate(condition_mean = r_Language) %>%
    ggplot(aes(y = Language, x = condition_mean, fill = after_stat(x < 0))) +
    stat_halfeye(alpha=0.75, point_interval=median_qi, .width=c(.50, .95)) +
    scale_fill_manual(values = c("lightsalmon", "skyblue"), name=NULL, labels=c("≥ 0", "< 0")) +
    geom_vline(xintercept = 0, linetype = "dashed") +
    ggtitle("1 | Language") + xlab("Intercept") + ylab(NULL) + xlim(-2,2) +
    theme(legend.position="bottom") +
    NULL,

  nrow=1);

Figure 28. Posterior estimaes of the random effects showing the median (dot), 50% (thick line) and 95% (thin lines) quantiles. Also showing 0 as the vertical dotted line and negative (blue) vs positive (red) values. Figure generated using R version 4.3.3 (2024-02-29)

Interpretation: The overall probability of a “perfect” match in the field experiment when adding all the relevant [r]-related predictors for the L1(s) and L2(s) known to the participants, while controlling for Language, is higher than the conservative chance level of 50%, being estimated as about 100% but with a much wider 95%CI of about [54%, 100.0%] which now includes 50%. No predictors seem to make any difference.

5 Both experiments together

Given that Order does not seem to formally make a difference in the online experiment, we can, in principle, combine the online and the field experiments in a joint analysis by ignoring Order as a predictor (but keeping all the observations).

webfield <- rbind(web[,names(field)], field); 
webfield$dataset <- factor(c(rep("web", nrow(web)), rep("field", nrow(field))), levels=c("web", "field"));
webfield$Sex <- toupper(webfield$Sex); # make sex uniform

5.1 What predictors to include?

We use the same ideas as for the web experiment above, except that now Order does not exist anymore.

webfield %>% count(Match) %>% mutate(prop = paste0(round(100* n / sum(n),1),"%")) 

webfield %>% count(r_l_distinction_L1) %>% mutate(prop = paste0(round(100* n / sum(n),1),"%")) 

webfield %>% count(trill_real_L1) %>% mutate(prop = paste0(round(100* n / sum(n),1),"%")) 

webfield %>% count(trill_occ_L1) %>% mutate(prop = paste0(round(100* n / sum(n),1),"%")) 

webfield %>% count(r_l_distinction_L2) %>% mutate(prop = paste0(round(100* n / sum(n),1),"%")) 

webfield %>% count(trill_real_L2) %>% mutate(prop = paste0(round(100* n / sum(n),1),"%")) 

webfield %>% count(trill_occ_L2) %>% mutate(prop = paste0(round(100* n / sum(n),1),"%")) 

It can be seen that, in our data, trill_occ_L1 and trill_occ_L2 are (almost) at 100% “yes” (for the non-missing data cases), which means that we will not include them in any further analyses.

5.2 What random structure to use?

unique(webfield[,c("Name", "Family", "Autotyp_Area")]) %>% arrange(Autotyp_Area, Family, Name)

r/l distinction in L1 (r_l_distinction_L1):

table(webfield$r_l_distinction_L1, webfield$Language)

##           
##             AL  AM  BE  BR  CN  DE  DK  EE  EN  ES  FA  FI  FR  GE  GR  HU  IT
##   same       0   0   0   0  46   0   0   0   0   0   0   0   0   0   0   0   0
##   distinct  10  20  55  20   0 104  18  43  39  36  21  18  57  15  42  34  52
##           
##             JP  KR  PA  PL  PT  RO  RU  SE  SR  TH  TR  VA  ZU
##   same      55  22   8   0   0   0   0   0   0   0   0   0  20
##   distinct   0   0   0  53  61  31  47  21  13  20  37  12   0

table(webfield$r_l_distinction_L1, webfield$Family)

##           
##            Afro-Asiatic Arawakan Austronesian Benue-Congo Finno-Ugric  IE
##   same                0        8            0          20           0   0
##   distinct           20        0           12           0          95 680
##           
##            Japanese Kartvelian Korean Sino-Tibetan Tai-Kadai Turkic
##   same           55          0     22           46         0      0
##   distinct        0         15      0            0        20     37

table(webfield$r_l_distinction_L1, webfield$Autotyp_Area)

##           
##            Europe Greater Mesopotamia Inner Asia N Africa N Coast Asia
##   same          0                   0          0        0           77
##   distinct    626                  93        108       20            0
##           
##            NE South America Oceania S Africa Southeast Asia
##   same                    8       0       20             46
##   distinct                0      12        0             20

This is perfectly uniform within the levels of each of the three factors, so random slopes are arguably not justified.

r/l distinction in L2 (r_l_distinction_L2):

table(webfield$r_l_distinction_L2, webfield$Language)

##           
##            AL AM BE BR CN DE DK EE EN ES FA FI FR GE GR HU IT JP KR PA PL PT RO
##   same      0  0  1  0  0  0  0  0  1  0  0  0  0  0  0  0  0  1  0  0  0  0  0
##   distinct  9 20 18 20 30 98 16 43 28 34 18 18 46 15 42 26 50 28 21  8 52 42 30
##           
##            RU SE SR TH TR VA ZU
##   same      0  0  0  0  0  0  0
##   distinct 43 19  1 18 29 12 18

table(webfield$r_l_distinction_L2, webfield$Family)

##           
##            Afro-Asiatic Arawakan Austronesian Benue-Congo Finno-Ugric  IE
##   same                0        0            0           0           0   2
##   distinct           20        8           12          18          87 566
##           
##            Japanese Kartvelian Korean Sino-Tibetan Tai-Kadai Turkic
##   same            1          0      0            0         0      0
##   distinct       28         15     21           30        18     29

table(webfield$r_l_distinction_L2, webfield$Autotyp_Area)

##           
##            Europe Greater Mesopotamia Inner Asia N Africa N Coast Asia
##   same          2                   0          0        0            1
##   distinct    511                  82        104       20           49
##           
##            NE South America Oceania S Africa Southeast Asia
##   same                    0       0        0              0
##   distinct                8      12       18             48

There is almost uniform, so random slopes are arguably not justified.

[r] is main allophone in L1 (trill_real_L1):

table(webfield$trill_real_L1, webfield$Language)

##      
##        AL  AM  BE  BR  CN  DE  DK  EE  EN  ES  FA  FI  FR  GE  GR  HU  IT  JP
##   no    0   0  55   0  46 101  18   0  38   0   0   0  56  14  42   0   0  55
##   yes  10  20   0  20   0   3   0  43   1  36  21  18   1   1   0  34  52   0
##      
##        KR  PA  PL  PT  RO  RU  SE  SR  TH  TR  VA  ZU
##   no   22   8   0  61   0   0  20  13  20  37   0  20
##   yes   0   0  53   0  31  47   1   0   0   0  12   0

table(webfield$trill_real_L1, webfield$Family)

##      
##       Afro-Asiatic Arawakan Austronesian Benue-Congo Finno-Ugric  IE Japanese
##   no             0        8            0          20           0 404       55
##   yes           20        0           12           0          95 276        0
##      
##       Kartvelian Korean Sino-Tibetan Tai-Kadai Turkic
##   no          14     22           46        20     37
##   yes          1      0            0         0      0

table(webfield$trill_real_L1, webfield$Autotyp_Area)

##      
##       Europe Greater Mesopotamia Inner Asia N Africa N Coast Asia
##   no     404                  51          0        0           77
##   yes    222                  42        108       20            0
##      
##       NE South America Oceania S Africa Southeast Asia
##   no                 8       0       20             66
##   yes                0      12        0              0

This is almost uniform within Language, but there is variation for the IE level of Family, and between 2 levels of Area (Europeand Greater Mesopotamia), but this is not enough to justify random slopes.

[r] is main allophone in L2 (trill_real_L2):

table(webfield$trill_real_L2, webfield$Language)

##      
##       AL AM BE BR CN DE DK EE EN ES FA FI FR GE GR HU IT JP KR PA PL PT RO RU
##   no   5  1  9  0 29 48 13  5 16 22 16  8 28  3 28 18 30 26 21  8 34 19 21 32
##   yes  4 19 10 20  1 50  3 38 13 12  2 10 18 12 14  8 20  3  0  0 18 23  9 11
##      
##       SE SR TH TR VA ZU
##   no  11  1 18 28  0 16
##   yes  8  0  0  1 12  2

table(webfield$trill_real_L2, webfield$Family)

##      
##       Afro-Asiatic Arawakan Austronesian Benue-Congo Finno-Ugric  IE Japanese
##   no             0        8            0          16          31 334       26
##   yes           20        0           12           2          56 234        3
##      
##       Kartvelian Korean Sino-Tibetan Tai-Kadai Turkic
##   no           3     21           29        18     28
##   yes         12      0            1         0      1

table(webfield$trill_real_L2, webfield$Autotyp_Area)

##      
##       Europe Greater Mesopotamia Inner Asia N Africa N Coast Asia
##   no     303                  48         45        0           47
##   yes    210                  34         59       20            3
##      
##       NE South America Oceania S Africa Southeast Asia
##   no                 8       0       16             47
##   yes                0      12        2              1

Here there is enough variation within the levels for all three factors, justifying the inclusion of random slopes.

Given these, we did the following:

• we do not include random slopes for the r/l distinction in L1 (r_l_distinction_L1) or L2 (r_l_distinction_L2), nor for the presence of [r] in L1 (trill_real_L1), but
• we do include random slopes for presence of [r] in L2 (trill_real_L2) by Language, Family and Area.

5.3 Model 1: what is the overall probability of a match?

if( !file.exists("./cached_results/b1_res_webfield.RData") )
{
  # prior predictive checks:
  # what priors we can set:
  get_prior(Match ~ 1 + dataset + 
              (1 | Language) +
              (1 | Family) +
              (1 | Autotyp_Area),
            data=webfield,
            family=bernoulli(link='logit'));
  b_priors <- brms::brm(Match ~ 1 + dataset + 
                          (1 | Language) +
                          (1 | Family) +
                          (1 | Autotyp_Area),
                        data=webfield,
                        family=bernoulli(link='logit'),
                        prior=c(brms::set_prior("student_t(5, 0, 2.5)", class="b")),
                        sample_prior='only',  # needed for prior predictive checks
                        seed=998, cores=brms_ncores, iter=10000, warmup=4000, thin=2, control=list(adapt_delta=0.999, max_treedepth=13));
  grid.arrange(pp_check(b_priors, type='stat', sta ='min') + xlab('p(match)') + ggtitle('Prior predictive distribution of minimum values'),
               pp_check(b_priors, type='stat', stat='mean') + xlab('p(match)') + ggtitle('Prior predictive distribution of means'),
               pp_check(b_priors, type='stat', stat='max') + xlab('p(match)') + ggtitle('Prior predictive distribution of maximum values'),
               ncol=1); # seems fine but our value is a bit extreme
  
  # the model:
  b1 <- brm(Match ~ 1 + dataset + 
              (1 | Language) +
              (1 | Family) +
              (1 | Autotyp_Area),
            data=webfield,
            family=bernoulli(link='logit'),
            prior=c(brms::set_prior("student_t(5, 0, 2.5)", class="b")),
            save_pars=save_pars(all=TRUE), # needed for Bayes factors
            sample_prior=TRUE,  # needed for hypotheses tests
            seed=998, cores=brms_ncores, iter=10000, warmup=4000, thin=2, control=list(adapt_delta=0.999999, max_treedepth=13));
  summary(b1); mcmc_plot(b1, type="trace"); mcmc_plot(b1); # very decent
  mcmc_plot(b1, type="trace");
  bayestestR::hdi(b1, ci=0.95);
  hypothesis(b1, "datasetfield = 0"); 
  # posterior predictive checks
  pp_check(b1, ndraws=100) + xlab('p(match)') + ggtitle('Posterior predictive density overlay');
  grid.arrange(pp_check(b1, type='stat', sta ='min') + xlab('p(match)') + ggtitle('Posterior predictive distribution of minimum values'),
               pp_check(b1, type='stat', stat='mean') + xlab('p(match)') + ggtitle('Posterior predictive distribution of means'),
               pp_check(b1, type='stat', stat='max') + xlab('p(match)') + ggtitle('Posterior predictive distribution of maximum values'),
               ncol=1); # far from ideal but seems ok...
  
  # save the model:
  b1_res_webfield <- list("b_priors"=b_priors, "b1"=b1);
  saveRDS(b1_res_webfield, "./cached_results/b1_res_webfield.RData", compress="xz");
} else
{
  b1_res_webfield <- readRDS("./cached_results/b1_res_webfield.RData");
}

grid.arrange(pp_check(b1_res_webfield$b_priors, type='stat', sta ='min') + xlab('p(match)') + ggtitle('Prior predictive distribution of minimum values'),
             pp_check(b1_res_webfield$b_priors, type='stat', stat='mean') + xlab('p(match)') + ggtitle('Prior predictive distribution of means'),
             pp_check(b1_res_webfield$b_priors, type='stat', stat='max') + xlab('p(match)') + ggtitle('Prior predictive distribution of maximum values'),
             ncol=1); # seems fine but our value is a bit extreme

## Using all posterior draws for ppc type 'stat' by default.
## Using all posterior draws for ppc type 'stat' by default.
## Using all posterior draws for ppc type 'stat' by default.

## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.

Figure 29. Prior predictive checks. Figure generated using R version 4.3.3 (2024-02-29)

and it can be seen that they are quite ok, even if the observed p(match) is a bit extreme, but still within what the prior distribution covers.

Convergence The model converges well (high ESS and \(\hat{R} = 1.0\)):

summary(b1_res_webfield$b1);

##  Family: bernoulli 
##   Links: mu = logit 
## Formula: Match ~ 1 + dataset + (1 | Language) + (1 | Family) + (1 | Autotyp_Area) 
##    Data: webfield (Number of observations: 1030) 
##   Draws: 4 chains, each with iter = 10000; warmup = 4000; thin = 2;
##          total post-warmup draws = 12000
## 
## Multilevel Hyperparameters:
## ~Autotyp_Area (Number of levels: 9) 
##               Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## sd(Intercept)     0.40      0.36     0.02     1.31 1.00     6971     8935
## 
## ~Family (Number of levels: 12) 
##               Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## sd(Intercept)     0.34      0.29     0.01     1.05 1.00     7936     9857
## 
## ~Language (Number of levels: 30) 
##               Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## sd(Intercept)     0.57      0.19     0.23     0.98 1.00     6446     7944
## 
## Regression Coefficients:
##              Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## Intercept        2.01      0.33     1.39     2.71 1.00    10336    10072
## datasetfield     1.64      0.68     0.39     3.06 1.00    10442    10343
## 
## Draws were sampled using sampling(NUTS). For each parameter, Bulk_ESS
## and Tail_ESS are effective sample size measures, and Rhat is the potential
## scale reduction factor on split chains (at convergence, Rhat = 1).

mcmc_plot(b1_res_webfield$b1, type="trace");

## No divergences to plot.

Figure 30. Trace of the fitting process. Figure generated using R version 4.3.3 (2024-02-29)

Posterior predictive checks seems to be acceptable:

pp_check(b1_res_webfield$b1, ndraws=500) + xlab('p(match)') + ggtitle('Posterior predictive density overlay');

Figure 31. Posterior predictive density overlay for n=500 draws. Figure generated using R version 4.3.3 (2024-02-29)

grid.arrange(pp_check(b1_res_webfield$b1, type='stat', sta ='min') + xlab('p(match)') + ggtitle('Posterior predictive distribution of minimum values') + xlim(0,1),
             pp_check(b1_res_webfield$b1, type='stat', stat='mean') + xlab('p(match)') + ggtitle('Posterior predictive distribution of means') + xlim(0,1),
             pp_check(b1_res_webfield$b1, type='stat', stat='max') + xlab('p(match)') + ggtitle('Posterior predictive distribution of maximum values') + xlim(0,1),
             ncol=1); 

## Using all posterior draws for ppc type 'stat' by default.
## Using all posterior draws for ppc type 'stat' by default.
## Using all posterior draws for ppc type 'stat' by default.

## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.

Figure 32. Posterior predictive distribution. Figure generated using R version 4.3.3 (2024-02-29)

Fixed effects:

mcmc_plot(b1_res_webfield$b1);

Figure 33. Fixed effects’ posterior distributions. Figure generated using R version 4.3.3 (2024-02-29)

plot_model(b1_res_webfield$b1, type="emm", terms=c("dataset"));

Figure 34. Fixed effects. Figure generated using R version 4.3.3 (2024-02-29)

Intercept The intercept is 2.01 95%CI [1.39, 2.71] on the log odds ratio (LOR) scale, which is clearly » 0; formal hypothesis testing against 0 is p(α>0)=1.00*; this translates into a probability of a match p(match)=88.2% 95%CI [80.0%, 93.8%] » 50%. (N.B. this is for dataset == “web”; please below for details.)

Dataset (experiment) The slope of dataset is 1.64 95%CI [0.39, 3.06] on the log odds ratio (LOR) scale, which is clearly > 0; formal hypothesis testing against 0 is p(β>0)=0.99*. In terms of probabilities, the field experiment (dataset == “field”) results in an increase by 9.3% in the probability of a match from the 88.2% [80.0%, 93.8%] in the web experiment, to 97.5% [91.5%, 99.4%] in the field experiment.

Random effects:

#get_variables(b1_res_webfield$b1);
grid.arrange(
  # 1 | Language:
  b1_res_webfield$b1 %>%
    spread_draws(b_Intercept, r_Language[Language, Order]) %>%
    filter(Order == "Intercept") %>%
    mutate(condition_mean = r_Language) %>%
    ggplot(aes(y = Language, x = condition_mean, fill = after_stat(x < 0))) +
    stat_halfeye(alpha=0.75, point_interval=median_qi, .width=c(.50, .95)) +
    scale_fill_manual(values = c("lightsalmon", "skyblue"), name=NULL, labels=c("≥ 0", "< 0")) +
    geom_vline(xintercept = 0, linetype = "dashed") +
    ggtitle("1 | Language") + xlab("Intercept") + ylab(NULL) + xlim(-2,2) +
    theme(legend.position="bottom") +
    NULL,
  
  # 1 | Family:
  b1_res_webfield$b1 %>%
    spread_draws(b_Intercept, r_Family[Family, Order]) %>%
    filter(Order == "Intercept") %>%
    mutate(condition_mean = r_Family) %>%
    ggplot(aes(y = Family, x = condition_mean, fill = after_stat(x < 0))) +
    stat_halfeye(alpha=0.75, point_interval=median_qi, .width=c(.50, .95)) +
    scale_fill_manual(values = c("lightsalmon", "skyblue"), name=NULL, labels=c("≥ 0", "< 0")) +
    geom_vline(xintercept = 0, linetype = "dashed") +
    ggtitle("1 | Family") + xlab("Intercept") + ylab(NULL) + xlim(-2,2) +
    theme(legend.position="bottom") +
    NULL,
  
  # 1 | Autotyp_Area:
  b1_res_webfield$b1 %>%
    spread_draws(b_Intercept, r_Autotyp_Area[Autotyp_Area, Order]) %>%
    filter(Order == "Intercept") %>%
    mutate(condition_mean = r_Autotyp_Area) %>%
    ggplot(aes(y = Autotyp_Area, x = condition_mean, fill = after_stat(x < 0))) +
    stat_halfeye(alpha=0.75, point_interval=median_qi, .width=c(.50, .95)) +
    scale_fill_manual(values = c("lightsalmon", "skyblue"), name=NULL, labels=c("≥ 0", "< 0")) +
    geom_vline(xintercept = 0, linetype = "dashed") +
    ggtitle("1 | Area") + xlab("Intercept") + ylab(NULL) + xlim(-2,2) +
    theme(legend.position="bottom") +
    NULL,

  nrow=1);

Figure 35. Posterior estimaes of the random effects showing the median (dot), 50% (thick line) and 95% (thin lines) quantiles. Also showing 0 as the vertical dotted line and negative (blue) vs positive (red) values. Figure generated using R version 4.3.3 (2024-02-29)

Interpretation: The overall probability of a “perfect” match, while controlling for Language, Family and Area, is much higher than the conservative change level of 50%, being estimated at more than 88% [80%, 94%] in the web experiment, and more than 97% [92%, 99%] in the field experiment, suggesting that there is a very strong tendency across our participants to associate [r] with the wavy line and [l] with the straight line.

5.4 Model 2: do the characteristics of the languages matter?

if( !file.exists("./cached_results/b2_res_webfield.RData") )
{
  # the model:
  b2 <- brm(Match ~ 1 + dataset + # dataset
              r_l_distinction_L1 + r_l_distinction_L2 + trill_real_L1 + trill_real_L2 + # main effects
              dataset:r_l_distinction_L1 + dataset:r_l_distinction_L2 + dataset:trill_real_L1 + dataset:trill_real_L2 + # potential interactions with dataset
              (1 + trill_real_L2 | Language) +
              (1 + trill_real_L2 | Family) +
              (1 + trill_real_L2 | Autotyp_Area),
            data=webfield,
            family=bernoulli(link='logit'),
            prior=c(brms::set_prior("student_t(5, 0, 2.5)", class="b"), 
                    brms::set_prior("lkj(2)", class="cor")),
            save_pars=save_pars(all=TRUE), # needed for Bayes factors
            sample_prior=TRUE,  # needed for hypotheses tests
            seed=998, cores=brms_ncores, iter=10000, warmup=4000, thin=2, control=list(adapt_delta=0.999, max_treedepth=13));
  summary(b2); mcmc_plot(b2, type="trace"); mcmc_plot(b2); # very decent
  mcmc_plot(b2, type="trace");
  bayestestR::hdi(b2, ci=0.95);
  hypothesis(b2, c("datasetfield = 0", 
                   "r_l_distinction_L1distinct = 0", "r_l_distinction_L2distinct = 0", "trill_real_L1yes = 0", "trill_real_L2yes = 0",
                   "datasetfield:r_l_distinction_L1distinct = 0", "datasetfield:r_l_distinction_L2distinct = 0", "datasetfield:trill_real_L1yes = 0", "datasetfield:trill_real_L2yes = 0")); 
  # posterior predictive checks
  pp_check(b2, ndraws=100) + xlab('p(match)') + ggtitle('Posterior predictive density overlay');
  grid.arrange(pp_check(b2, type='stat', sta ='min') + xlab('p(match)') + ggtitle('Posterior predictive distribution of minimum values') + xlim(0,1),
               pp_check(b2, type='stat', stat='mean') + xlab('p(match)') + ggtitle('Posterior predictive distribution of means') + xlim(0,1),
               pp_check(b2, type='stat', stat='max') + xlab('p(match)') + ggtitle('Posterior predictive distribution of maximum values') + xlim(0,1),
               ncol=1); # seems perfect
  
  # the interactions of dataset do not seem to matter at all, so let's remove them:
  b3 <- brm(Match ~ 1 + dataset + # dataset
              r_l_distinction_L1 + r_l_distinction_L2 + trill_real_L1 + trill_real_L2 + # main effects
              (1 + trill_real_L2 | Language) +
              (1 + trill_real_L2 | Family) +
              (1 + trill_real_L2 | Autotyp_Area),
            data=webfield,
            family=bernoulli(link='logit'),
            prior=c(brms::set_prior("student_t(5, 0, 2.5)", class="b"), 
                    brms::set_prior("lkj(2)", class="cor")),
            save_pars=save_pars(all=TRUE), # needed for Bayes factors
            sample_prior=TRUE,  # needed for hypotheses tests
            seed=998, cores=brms_ncores, iter=10000, warmup=4000, thin=2, control=list(adapt_delta=0.999, max_treedepth=13));
  summary(b3); mcmc_plot(b3, type="trace"); mcmc_plot(b3); # very decent
  mcmc_plot(b3, type="trace");
  bayestestR::hdi(b3, ci=0.95);
  hypothesis(b3, c("datasetfield = 0", "r_l_distinction_L1distinct = 0", "r_l_distinction_L2distinct = 0", "trill_real_L1yes = 0", "trill_real_L2yes = 0")); 
  # posterior predictive checks
  pp_check(b3, ndraws=100) + xlab('p(match)') + ggtitle('Posterior predictive density overlay');
  grid.arrange(pp_check(b3, type='stat', sta ='min') + xlab('p(match)') + ggtitle('Posterior predictive distribution of minimum values') + xlim(0,1),
               pp_check(b3, type='stat', stat='mean') + xlab('p(match)') + ggtitle('Posterior predictive distribution of means') + xlim(0,1),
               pp_check(b3, type='stat', stat='max') + xlab('p(match)') + ggtitle('Posterior predictive distribution of maximum values') + xlim(0,1),
               ncol=1); # seems perfect
  
  # save the model:
  b2_res_webfield <- list("b2"=b2, "b3"=b3);
  saveRDS(b2_res_webfield, "./cached_results/b2_res_webfield.RData", compress="xz");
} else
{
  b2_res_webfield <- readRDS("./cached_results/b2_res_webfield.RData");
}

We fist included the interactions between dataset and the [r]-related predictors, but they are not relevant (it is interesting to note that even in this model, trill_real_L1 has a significant effect):

bayestestR::hdi(b2_res_webfield$b2, ci=0.95);

hypothesis(b2_res_webfield$b2, c("datasetfield = 0", 
                                 "r_l_distinction_L1distinct = 0", "r_l_distinction_L2distinct = 0", "trill_real_L1yes = 0", "trill_real_L2yes = 0",
                                 "datasetfield:r_l_distinction_L1distinct = 0", "datasetfield:r_l_distinction_L2distinct = 0", "datasetfield:trill_real_L1yes = 0", "datasetfield:trill_real_L2yes = 0"));

## Hypothesis Tests for class b:
##                 Hypothesis Estimate Est.Error CI.Lower CI.Upper Evid.Ratio
## 1       (datasetfield) = 0     2.49      2.95    -2.31     9.31       0.80
## 2 (r_l_distinction_... = 0     0.49      0.93    -1.47     2.25       2.55
## 3 (r_l_distinction_... = 0    -1.29      2.59    -7.11     2.59       1.26
## 4   (trill_real_L1yes) = 0    -0.99      0.40    -1.81    -0.24       0.27
## 5   (trill_real_L2yes) = 0     0.38      0.93    -1.27     2.47       3.52
## 6 (datasetfield:r_l... = 0     1.50      2.73    -3.17     7.78       1.05
## 7 (datasetfield:r_l... = 0     2.24      2.86    -2.46     8.74       0.88
## 8 (datasetfield:tri... = 0     0.80      2.78    -4.31     6.84       1.09
## 9 (datasetfield:tri... = 0     0.95      2.86    -4.09     6.88       1.12
##   Post.Prob Star
## 1      0.44     
## 2      0.72     
## 3      0.56     
## 4      0.21    *
## 5      0.78     
## 6      0.51     
## 7      0.47     
## 8      0.52     
## 9      0.53     
## ---
## 'CI': 90%-CI for one-sided and 95%-CI for two-sided hypotheses.
## '*': For one-sided hypotheses, the posterior probability exceeds 95%;
## for two-sided hypotheses, the value tested against lies outside the 95%-CI.
## Posterior probabilities of point hypotheses assume equal prior probabilities.

so we removed them, keeping only the main effects of the predictors (including dataset).

Priors are as before.

Convergence The model converges well (high ESS and \(\hat{R} = 1.0\)):

summary(b2_res_webfield$b3);

##  Family: bernoulli 
##   Links: mu = logit 
## Formula: Match ~ 1 + dataset + r_l_distinction_L1 + r_l_distinction_L2 + trill_real_L1 + trill_real_L2 + (1 + trill_real_L2 | Language) + (1 + trill_real_L2 | Family) + (1 + trill_real_L2 | Autotyp_Area) 
##    Data: webfield (Number of observations: 855) 
##   Draws: 4 chains, each with iter = 10000; warmup = 4000; thin = 2;
##          total post-warmup draws = 12000
## 
## Multilevel Hyperparameters:
## ~Autotyp_Area (Number of levels: 9) 
##                                 Estimate Est.Error l-95% CI u-95% CI Rhat
## sd(Intercept)                       0.57      0.50     0.02     1.89 1.00
## sd(trill_real_L2yes)                1.07      0.93     0.04     3.48 1.00
## cor(Intercept,trill_real_L2yes)    -0.03      0.45    -0.83     0.81 1.00
##                                 Bulk_ESS Tail_ESS
## sd(Intercept)                       6674     9320
## sd(trill_real_L2yes)                7321     8492
## cor(Intercept,trill_real_L2yes)     9759     9994
## 
## ~Family (Number of levels: 12) 
##                                 Estimate Est.Error l-95% CI u-95% CI Rhat
## sd(Intercept)                       0.93      0.56     0.06     2.19 1.00
## sd(trill_real_L2yes)                0.95      0.89     0.03     3.22 1.00
## cor(Intercept,trill_real_L2yes)    -0.16      0.45    -0.88     0.74 1.00
##                                 Bulk_ESS Tail_ESS
## sd(Intercept)                       5023     6050
## sd(trill_real_L2yes)                7748     9334
## cor(Intercept,trill_real_L2yes)    10006    10648
## 
## ~Language (Number of levels: 30) 
##                                 Estimate Est.Error l-95% CI u-95% CI Rhat
## sd(Intercept)                       0.38      0.27     0.02     0.99 1.00
## sd(trill_real_L2yes)                0.87      0.48     0.08     1.94 1.00
## cor(Intercept,trill_real_L2yes)    -0.10      0.44    -0.84     0.75 1.00
##                                 Bulk_ESS Tail_ESS
## sd(Intercept)                       5731     7352
## sd(trill_real_L2yes)                5253     6479
## cor(Intercept,trill_real_L2yes)     6527     9138
## 
## Regression Coefficients:
##                            Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS
## Intercept                      3.27      2.57    -0.82     9.14 1.00     8579
## datasetfield                   4.22      2.58     0.92    10.16 1.00     8095
## r_l_distinction_L1distinct     0.54      0.93    -1.45     2.33 1.00     8002
## r_l_distinction_L2distinct    -1.21      2.47    -6.85     2.59 1.00     8657
## trill_real_L1yes              -0.98      0.41    -1.81    -0.23 1.00     8916
## trill_real_L2yes               0.47      0.95    -1.12     2.70 1.00     8081
##                            Tail_ESS
## Intercept                      8945
## datasetfield                   6470
## r_l_distinction_L1distinct     8899
## r_l_distinction_L2distinct     8438
## trill_real_L1yes               9618
## trill_real_L2yes               8436
## 
## Draws were sampled using sampling(NUTS). For each parameter, Bulk_ESS
## and Tail_ESS are effective sample size measures, and Rhat is the potential
## scale reduction factor on split chains (at convergence, Rhat = 1).

mcmc_plot(b2_res_webfield$b3, type="trace");

## No divergences to plot.

Figure 36. Trace of the fitting process. Figure generated using R version 4.3.3 (2024-02-29)

Posterior predictive checks seems to be perfectly fine:

pp_check(b2_res_webfield$b3, ndraws=500) + xlab('p(match)') + ggtitle('Posterior predictive density overlay');

Figure 37. Posterior predictive density overlay for n=500 draws. Figure generated using R version 4.3.3 (2024-02-29)

grid.arrange(pp_check(b2_res_webfield$b3, type='stat', sta ='min') + xlab('p(match)') + ggtitle('Posterior predictive distribution of minimum values') + xlim(0,1),
             pp_check(b2_res_webfield$b3, type='stat', stat='mean') + xlab('p(match)') + ggtitle('Posterior predictive distribution of means') + xlim(0,1),
             pp_check(b2_res_webfield$b3, type='stat', stat='max') + xlab('p(match)') + ggtitle('Posterior predictive distribution of maximum values') + xlim(0,1),
             ncol=1); 

## Using all posterior draws for ppc type 'stat' by default.
## Using all posterior draws for ppc type 'stat' by default.
## Using all posterior draws for ppc type 'stat' by default.

## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.

Figure 38. Posterior predictive distribution (the x-axis scale was forced to [0,1]). Figure generated using R version 4.3.3 (2024-02-29)

Fixed effects:

mcmc_plot(b2_res_webfield$b3);

Figure 39. Fixed effects’ posterior distributions. Figure generated using R version 4.3.3 (2024-02-29)

grid.arrange(plot_model(b2_res_webfield$b3, type="emm", terms=c("dataset")),
             plot_model(b2_res_webfield$b3, type="emm", terms=c("r_l_distinction_L1")),
             plot_model(b2_res_webfield$b3, type="emm", terms=c("r_l_distinction_L2")),
             plot_model(b2_res_webfield$b3, type="emm", terms=c("trill_real_L1")),
             plot_model(b2_res_webfield$b3, type="emm", terms=c("trill_real_L2")),
             ncol=2);

Figure 40. Fixed effects. Figure generated using R version 4.3.3 (2024-02-29)

Intercept The intercept is 3.27 95%CI [-0.82, 9.14] on the log odds ratio (LOR) scale, which is > 0 but has very large uncertainties so that the 95%CI includes 0; formal hypothesis testing against 0 is p(α>0)=0.93; this translates into a probability of a match p(match)=96.4% 95%CI [30.6%, 100.0%] which includes 50% in the 95%CI.

Dataset (experiment) The slope of dataset is 4.22 95%CI [0.92, 10.16] on the log odds ratio (LOR) scale, which is clearly > 0; formal hypothesis testing against 0 is p(β>0)=1.00*. In terms of probabilities, the field experiment (dataset == “field”) results in an increase by 3.6% in the probability of a match from the 96.4% [30.6%, 100.0%] in the web experiment, to 99.9% [87.8%, 100.0%] in the field experiment.

Are [r] and [l] distinct in the L1(s) spoken by the participant? The slope of r_l_distinction_L1 is 0.54 95%CI [-1.45, 2.33] on the log odds ratio (LOR) scale, which includes 0; formal hypothesis testing against 0 is p(β>0)=0.75.

Are [r] and [l] distinct in the L2(s) spoken by the participant? The slope of r_l_distinction_L2 is -1.21 95%CI [-6.85, 2.59] on the log odds ratio (LOR) scale, which includes 0; formal hypothesis testing against 0 is p(β<0)=0.67.

Is [r] the main allophone in the L1(s) spoken by the participant? The slope of trill_real_L1 is -0.98 95%CI [-1.81, -0.23] on the log odds ratio (LOR) scale, which is clearly < 0; formal hypothesis testing against 0 is p(β<0)=0.99*. In terms of probabilities, having [r] as the main allophone (trill_real_L1 == “yes”) results in a decrease by 5.6% [0.003%, 23.9%] in the probability of a match from the 96.4% [30.6%, 100.0%] when [r] is not the primary allophone to 90.8% [13.7%, 100.0%] when it is.

Is [r] the main allophone in the L2(s) spoken by the participant? The slope of trill_real_L2 is 0.47 95%CI [-1.12, 2.70] on the log odds ratio (LOR) scale, which includes 0; formal hypothesis testing against 0 with p(β>0)=0.68.

Random effects:

#get_variables(b2_res_webfield$b3);
grid.arrange(
  # 1 | Language:
  b2_res_webfield$b3 %>%
    spread_draws(b_Intercept, r_Language[Language, Order]) %>%
    filter(Order == "Intercept") %>%
    mutate(condition_mean = r_Language) %>%
    ggplot(aes(y = Language, x = condition_mean, fill = after_stat(x < 0))) +
    stat_halfeye(alpha=0.75, point_interval=median_qi, .width=c(.50, .95)) +
    scale_fill_manual(values = c("lightsalmon", "skyblue"), name=NULL, labels=c("≥ 0", "< 0")) +
    geom_vline(xintercept = 0, linetype = "dashed") +
    ggtitle("1 | Language") + xlab("Intercept") + ylab(NULL) + xlim(-2,2) +
    theme(legend.position="bottom") +
    NULL,
  
  # 1 | Family:
  b2_res_webfield$b3 %>%
    spread_draws(b_Intercept, r_Family[Family, Order]) %>%
    filter(Order == "Intercept") %>%
    mutate(condition_mean = r_Family) %>%
    ggplot(aes(y = Family, x = condition_mean, fill = after_stat(x < 0))) +
    stat_halfeye(alpha=0.75, point_interval=median_qi, .width=c(.50, .95)) +
    scale_fill_manual(values = c("lightsalmon", "skyblue"), name=NULL, labels=c("≥ 0", "< 0")) +
    geom_vline(xintercept = 0, linetype = "dashed") +
    ggtitle("1 | Family") + xlab("Intercept") + ylab(NULL) + xlim(-2,2) +
    theme(legend.position="bottom") +
    NULL,
  
  # 1 | Autotyp_Area:
  b2_res_webfield$b3 %>%
    spread_draws(b_Intercept, r_Autotyp_Area[Autotyp_Area, Order]) %>%
    filter(Order == "Intercept") %>%
    mutate(condition_mean = r_Autotyp_Area) %>%
    ggplot(aes(y = Autotyp_Area, x = condition_mean, fill = after_stat(x < 0))) +
    stat_halfeye(alpha=0.75, point_interval=median_qi, .width=c(.50, .95)) +
    scale_fill_manual(values = c("lightsalmon", "skyblue"), name=NULL, labels=c("≥ 0", "< 0")) +
    geom_vline(xintercept = 0, linetype = "dashed") +
    ggtitle("1 | Area") + xlab("Intercept") + ylab(NULL) + xlim(-2,2) +
    theme(legend.position="bottom") +
    NULL,

  # trill_real_L2 | Language:
  b2_res_webfield$b3 %>%
    spread_draws(b_Intercept, r_Language[Language, trill_real_L2]) %>%
    filter(trill_real_L2 == "trill_real_L2yes") %>%
    mutate(condition_mean = r_Language) %>%
    ggplot(aes(y = Language, x = condition_mean, fill = after_stat(x < 0))) +
    stat_halfeye(alpha=0.75, point_interval=median_qi, .width=c(.50, .95)) +
    scale_fill_manual(values = c("lightsalmon", "skyblue"), name=NULL, labels=c("≥ 0", "< 0")) +
    geom_vline(xintercept = 0, linetype = "dashed") +
    ggtitle("trill_real_L2 | Language") + xlab("Slope") + ylab(NULL) + xlim(-4, 4) +
    theme(legend.position="bottom") +
    NULL,
  
  # trill_real_L2 | Family:
  b2_res_webfield$b3 %>%
    spread_draws(b_Intercept, r_Family[Family, trill_real_L2]) %>%
    filter(trill_real_L2 == "trill_real_L2yes") %>%
    mutate(condition_mean = r_Family) %>%
    ggplot(aes(y = Family, x = condition_mean, fill = after_stat(x < 0))) +
    stat_halfeye(alpha=0.75, point_interval=median_qi, .width=c(.50, .95)) +
    scale_fill_manual(values = c("lightsalmon", "skyblue"), name=NULL, labels=c("≥ 0", "< 0")) +
    geom_vline(xintercept = 0, linetype = "dashed") +
    ggtitle("trill_real_L2 | Family") + xlab("Slope") + ylab(NULL) + xlim(-3, 3) + 
    theme(legend.position="bottom") +
    NULL,
  
  # trill_real_L2 | Autotyp_Area:
  b2_res_webfield$b3 %>%
    spread_draws(b_Intercept, r_Autotyp_Area[Autotyp_Area, trill_real_L2]) %>%
    filter(trill_real_L2 == "trill_real_L2yes") %>%
    mutate(condition_mean = r_Autotyp_Area) %>%
    ggplot(aes(y = Autotyp_Area, x = condition_mean, fill = after_stat(x < 0))) +
    stat_halfeye(alpha=0.75, point_interval=median_qi, .width=c(.50, .95)) +
    scale_fill_manual(values = c("lightsalmon", "skyblue"), name=NULL, labels=c("≥ 0", "< 0")) +
    geom_vline(xintercept = 0, linetype = "dashed") +
    ggtitle("trill_real_L2 | Area") + xlab("Slope") + ylab(NULL) + xlim(-3, 3) +  
    theme(legend.position="bottom") +
    NULL,
  
  nrow=2);

Figure 41. Posterior estimaes of the random effects showing the median (dot), 50% (thick line) and 95% (thin lines) quantiles. Also showing 0 as the vertical dotted line and negative (blue) vs positive (red) values. Figure generated using R version 4.3.3 (2024-02-29)

Interpretation: The overall probability of a “perfect” match in the online experiment when adding all the relevant [r]-related predictors for the L1(s) and L2(s) known to the participants, while controlling for Language, Family and Area, seem to still be higher than the conservative change level of 50%, being estimated as about 96% but with a much wider 95% credible interval of about [31%, 100.0%] which now includes 50%. If [r] is the main allophone in the L1(s) known by the participant or not clearly matters, in that those participants that have an L1 with [r] as the main allophone have a slightly smaller probability of a perfect match of about 91% [14%, 100%] than the others. Also, there is a clear difference between the two datasets, with the web participants having a lower probability of a match than the field participants by about 4%, to 100% [88%, 100%].

6 Plots

Here we build the various plots for the paper.

6.1 Web experiment

p <- gather_draws(b2_res_web$b2, `b_.*`, regex = TRUE) %>%
  # remove the b_ for plotting
  #mutate(.variable = gsub("b_","", .variable)) %>%
  ggplot(aes(x = .value, y = .variable, fill = after_stat(x < 0))) +
  geom_vline(xintercept=0.0, linetype="dotted", color="gray30") +
  stat_halfeye(alpha=0.75, point_interval=median_qi, .width=c(.50, .95)) +
  scale_fill_manual(values = c("lightsalmon", "skyblue"), name=NULL, labels=c("≥ 0", "< 0")) +
  xlim(-10, 10) +
  scale_y_discrete(labels=c("b_Intercept" = "Intercept α", 
                            "b_Order1" = "order ([r] first vs [l] first)",
                            "b_r_l_distinction_L1distinct" = "are [r] and [l] contrasting in L1(s)?",
                            "b_trill_real_L1yes" = "is [r] the main allophone in L1(s)?",
                            "b_r_l_distinction_L2distinct" = "are [r] and [l] contrasting in L2(s)?",
                            "b_trill_real_L2yes" = "is [r] the main allophone in L2(s)?")) +
  ylab(NULL) + xlab("Estimate") + 
  theme_timo + # Tweak cosmetics:
  theme(axis.title.x = element_text(margin = margin(t = 12, b = 0,
                                                    r = 0, l = 0),
                                    face = 'bold', size = 14),
        axis.text.x = element_text(face = 'bold', size = 12),
        axis.text.y = element_text(face = 'bold', size = 12))+
  NULL;
# Show and save:
p;

Figure 42. Posterior estimaes of the fixed effects showing the median (dot), 50% (thick line) and 95% (thin lines) quantiles. Also showing 0 as the vertical dotted line and negative (blue) vs positive (red) values. Figure generated using R version 4.3.3 (2024-02-29)

ggsave(plot = p, filename = "./plots/figure_model_fixedeffects_web.pdf", width = 8, height = 4);
ggsave(plot = p, filename = "./plots/figure_model_fixedeffects_web.jpg", width = 8, height = 4);
ggsave(plot = p, filename = "./plots/figure_model_fixedeffects_web.tif", width = 8, height = 4, compression="lzw", dpi=600);

# take2: generate an exhaustive new dataset:
# starting from the actual language coding:
languages_data <- read.csv("./data/languages_data.csv", sep=","); # read the original language info
newdata_web <- data.frame(Name = tolower(unique(web$Name)));
newdata_web <- merge(newdata_web, languages_data, by.x="Name", by.y="Languages", all.x=TRUE, all.y=FALSE);
# L1 coding:
names(newdata_web)[3:5] <- c("r_l_distinction_L1", "trill_real_L1", "trill_occ_L1");
newdata_web$r_l_distinction_L1 <- factor(c("same", "distinct")[newdata_web$r_l_distinction_L1 + 1], levels=c("same", "distinct"));
newdata_web$trill_real_L1 <- factor(c("no", "yes")[newdata_web$trill_real_L1 + 1], levels=c("no", "yes"));
newdata_web$trill_occ_L1 <- factor(c("no", "yes")[newdata_web$trill_occ_L1 + 1], levels=c("no", "yes"));
# family and area:
newdata_web <- merge(newdata_web, unique(web[,c("Language", "Name", "Family", "Autotyp_Area")]) %>% mutate(Name = tolower(Name)), by="Name", all.x=TRUE, all.y=FALSE);
# add it to the both newdata:
newdata_webfield <- newdata_web[  ,c("Name", "Language", "Family", "Autotyp_Area", "r_l_distinction_L1", "trill_real_L1", "trill_occ_L1")];

# fit a new model considering only the L1 info:
if( !file.exists("./cached_results/b_plotting_web.RData") )
{
  b_plotting_web <- brm(Match ~ 1 + 
                          r_l_distinction_L1 + trill_real_L1 + 
                          (1 | Language) +
                          (1 | Family) +
                          (1 | Autotyp_Area),
                        data=web,
                        family=bernoulli(link='logit'),
                        prior=c(brms::set_prior("student_t(5, 0, 2.5)", class="b")),
                        save_pars=save_pars(all=TRUE), # needed for Bayes factors
                        sample_prior=TRUE,  # needed for hypotheses tests
                        seed=998, cores=brms_ncores, iter=10000, warmup=4000, thin=2, control=list(adapt_delta=0.999, max_treedepth=13));
  summary(b_plotting_web); mcmc_plot(b_plotting_web, type="trace"); mcmc_plot(b_plotting_web); # very decent
  bayestestR::hdi(b_plotting_web, ci=0.95);

  # save the model:
  saveRDS(b_plotting_web, "./cached_results/b_plotting_web.RData", compress="xz");
} else
{
  b_plotting_web <- readRDS("./cached_results/b_plotting_web.RData");
}
fit <- fitted(b_plotting_web, newdata=newdata_web, re_formula=NULL, robust=TRUE);
colnames(fit) = c('fit', 'se', 'lwr', 'upr');
newdata_web <- cbind(newdata_web, fit);

# Order predictions by descriptive average:
newdata_web <- arrange(newdata_web, fit);
newdata_web <- mutate(newdata_web, Language = factor(Language, levels = newdata_web$Language));

# How many languages are over 0.5? (will be reported in paper)
sum(newdata_web$lwr > 0.5);

# Finally, add the averages to the plot:
newdata_web$avg <- web_avg[match(newdata_web$Language, web_avg$Language), ]$M;

# Match language names into there and order:
newdata_web$Language <- web[match(newdata_web$Language, web$Language), ]$Name;
newdata_web[newdata_web$Language == 'Chinese', ]$Language <- 'Mandarin Chinese';
newdata_web <- mutate(newdata_web, Language = factor(Language, levels = newdata_web$Language));

# Setup the plot:
# Aesthetics and geom:
plot_web <- newdata_web %>% 
  ggplot(aes(x = Language, col = r_l_distinction_L1, y = fit, ymin = lwr, ymax = upr, shape = trill_real_L1)) +
  #scale_shape_manual(values = c("no" = 1, "yes" = 3)) +
  geom_errorbar(aes(col = r_l_distinction_L1), linewidth = 1.0, width = 0.3) +
  geom_point(size = 5) +
  geom_hline(yintercept = 0.5, linetype = 2, linewidth = 1.5, col = 'grey') + 
  geom_point(aes(y = avg), col = 'black', shape = 23, size = 5, stroke = 1.0, alpha = 0.75)  +
  ylim(0.5,1.00) +
  labs(x = '', y = 'Proportion of congruent responses') + # Axis labels:
  #ggtitle('Posterior medians and descriptive averages of congruent responses\nby language and R/L contrast (color) and [r] as primary rhotic in L1 (shape)') +
  scale_color_manual(values = c("same"=colorBlindBlack8[2], "distinct"=colorBlindBlack8[3]), labels=c("same"="no [r]/[l] contrast", "distinct"="[r]/[l] contrast")) + # Tweak cosmetics:
  scale_shape(labels=c("yes"="[r] is main allophone", "no"="[r] is not main allophone")) + 
  theme_timo + 
  theme(axis.text.x = element_text(angle = 45, hjust = 1, face = 'bold', size = 10),
        axis.text.y = element_text(face = 'bold', size = 10),
        axis.title = element_text(face = 'bold', size = 12),
        axis.title.y = element_text(face = 'bold', size = 12, margin = margin(t = 0, r = 35, b = 0, l = 0)),
        plot.title = element_text(face = 'bold', size = 14, margin = margin(t = 0, r = 0, b = 30, l = 0)),
        legend.text = element_text(size = 10),
        legend.title = element_blank(),
        legend.position = c(0.95, 0.10),
        legend.justification = c('right', 'bottom'))

# Save:
ggsave(plot = plot_web, filename = "./plots/figure_model_languages_web.pdf", width = 12, height = 6);
ggsave(plot = plot_web, filename = "./plots/figure_model_languages_web.jpg", width = 12, height = 6);
ggsave(plot = plot_web, filename = "./plots/figure_model_languages_web.tif", width = 12, height = 6, compression="lzw", dpi=600);

plot_web;

Figure 43. The distribution of the proportion of matches per language. Please note that this was obtained by fitting a model that includes only the relevant L1 characteristics, Match ~ 1 + r_l_distinction_L1 + trill_real_L1 + (1 | Language) + (1 | Family) + (1 | Autotyp_Area). Figure generated using R version 4.3.3 (2024-02-29)

6.2 Field experiment

p <- gather_draws(b2_res_field$b2, `b_.*`, regex = TRUE) %>%
  # remove the b_ for plotting
  #mutate(.variable = gsub("b_","", .variable)) %>%
  ggplot(aes(x = .value, y = .variable, fill = after_stat(x < 0))) +
  geom_vline(xintercept=0.0, linetype="dotted", color="gray30") +
  stat_halfeye(alpha=0.75, point_interval=median_qi, .width=c(.50, .95)) +
  scale_fill_manual(values = c("lightsalmon", "skyblue"), name=NULL, labels=c("≥ 0", "< 0")) +
  xlim(-10, 20) +
  scale_y_discrete(labels=c("b_Intercept" = "Intercept α", 
                            "b_r_l_distinction_L1distinct" = "are [r] and [l] contrasting in L1(s)?",
                            "b_trill_real_L1yes" = "is [r] the main allophone in L1(s)?",
                            "b_r_l_distinction_L2distinct" = "are [r] and [l] contrasting in L2(s)?",
                            "b_trill_real_L2yes" = "is [r] the main allophone in L2(s)?",
                            "b_trill_occ_L1yes" = "is [r] an allophone in L1(s)?")) +
  ylab(NULL) + xlab("Estimate") + 
  theme_timo + # Tweak cosmetics:
  theme(axis.title.x = element_text(margin = margin(t = 12, b = 0,
                                                    r = 0, l = 0),
                                    face = 'bold', size = 14),
        axis.text.x = element_text(face = 'bold', size = 12),
        axis.text.y = element_text(face = 'bold', size = 12))+
  NULL;
# Show and save:
p;

Figure 44. Posterior estimaes of the fixed effects showing the median (dot), 50% (thick line) and 95% (thin lines) quantiles. Also showing 0 as the vertical dotted line and negative (blue) vs positive (red) values. Figure generated using R version 4.3.3 (2024-02-29)

ggsave(plot = p, filename = "./plots/full_model_field_coefficients.pdf", width = 8, height = 4);
ggsave(plot = p, filename = "./plots/full_model_field_coefficients.jpg", width = 8, height = 4);
ggsave(plot = p, filename = "./plots/full_model_field_coefficients.tif", width = 8, height = 4, compression="lzw", dpi=600);

# take2: generate an exhaustive new dataset:
# starting from the actual language coding:
languages_data <- read.csv("./data/languages_data.csv", sep=","); # read the original language info
newdata_field <- data.frame(Name = tolower(unique(field$Name)));
newdata_field$Name2 <- newdata_field$Name;
newdata_field$Name[ newdata_field$Name == "brazilian portuguese" ] <- "portuguese";
newdata_field$Name[ newdata_field$Name == "english uk" ] <- "english";
newdata_field$Name[ newdata_field$Name == "tashlhiyt berber" ] <- "berber";
newdata_field <- merge(newdata_field, languages_data, by.x="Name", by.y="Languages", all.x=TRUE, all.y=FALSE);
# manually fix some entries:
# L1 coding:
names(newdata_field)[4:6] <- c("r_l_distinction_L1", "trill_real_L1", "trill_occ_L1");
newdata_field$r_l_distinction_L1 <- factor(c("same", "distinct")[newdata_field$r_l_distinction_L1 + 1], levels=c("same", "distinct"));
newdata_field$trill_real_L1 <- factor(c("no", "yes")[newdata_field$trill_real_L1 + 1], levels=c("no", "yes"));
newdata_field$trill_occ_L1 <- factor(c("no", "yes")[newdata_field$trill_occ_L1 + 1], levels=c("no", "yes"));
# family and area:
newdata_field <- merge(newdata_field, unique(field[,c("Language", "Name", "Family", "Autotyp_Area")]) %>% mutate(Name = tolower(Name)), by.x="Name2", by.y="Name", all.x=TRUE, all.y=FALSE);
# add it to the both newdata:
newdata_webfield <- rbind(newdata_webfield,
                          newdata_field[,c("Name", "Language", "Family", "Autotyp_Area", "r_l_distinction_L1", "trill_real_L1", "trill_occ_L1")]);

# fit a new model considering only the L1 info:
if( !file.exists("./cached_results/b_plotting_field.RData") )
{
  b_plotting_field <- brm(Match ~ 1 + 
                          r_l_distinction_L1 + trill_real_L1 + trill_occ_L1 +
                          (1 | Language),
                        data=field,
                        family=bernoulli(link='logit'),
                        prior=c(brms::set_prior("student_t(5, 0, 2.5)", class="b")),
                        save_pars=save_pars(all=TRUE), # needed for Bayes factors
                        sample_prior=TRUE,  # needed for hypotheses tests
                        seed=998, cores=brms_ncores, iter=10000, warmup=4000, thin=2, control=list(adapt_delta=0.999, max_treedepth=13));
  summary(b_plotting_field); mcmc_plot(b_plotting_field, type="trace"); mcmc_plot(b_plotting_field); # very decent
  bayestestR::hdi(b_plotting_field, ci=0.95);
 
  # save the model:
  saveRDS(b_plotting_field, "./cached_results/b_plotting_field.RData", compress="xz");
} else
{
  b_plotting_field <- readRDS("./cached_results/b_plotting_field.RData");
}
fit <- fitted(b_plotting_field, newdata=newdata_field, re_formula=NULL, robust=TRUE);
colnames(fit) = c('fit', 'se', 'lwr', 'upr');
newdata_field <- cbind(newdata_field, fit);

# Order predictions by descriptive average:
newdata_field <- arrange(newdata_field, fit);
newdata_field <- mutate(newdata_field, Language = factor(Language, levels = newdata_field$Language));

# How many languages are over 0.5? (will be reported in paper)
sum(newdata_field$lwr > 0.5);

# Finally, add the averages to the plot:
newdata_field$avg <- field_avg[match(newdata_field$Language, field_avg$Language), ]$M;

# Match language names into there and order:
newdata_field$Language <- field[match(newdata_field$Language, field$Language), ]$Name;
newdata_field <- mutate(newdata_field, Language = factor(Language, levels = newdata_field$Language));

# Setup the plot:
# Aesthetics and geom:
plot_field <- newdata_field %>% 
  ggplot(aes(x = Language, col = r_l_distinction_L1, y = fit, ymin = lwr, ymax = upr, shape = trill_real_L1)) +
  #scale_shape_manual(values = c("no" = 1, "yes" = 3)) +
  geom_errorbar(aes(col = r_l_distinction_L1), linewidth = 1.0, width = 0.3) +
  geom_point(size = 5) +
  geom_hline(yintercept = 0.5, linetype = 2, linewidth = 1.5, col = 'grey') + 
  geom_point(aes(y = avg), col = 'black', shape = 23, size = 5, stroke = 1.0, alpha = 0.75)  +
  ylim(0.5,1.00) +
  labs(x = '', y = 'Proportion of congruent responses') + # Axis labels:
  #ggtitle('Posterior medians and descriptive averages of congruent responses\nby language and R/L contrast (color) and [r] as primary rhotic in L1 (shape)') +
  scale_color_manual(values = c("same"=colorBlindBlack8[2], "distinct"=colorBlindBlack8[3]), labels=c("same"="no [r]/[l] contrast", "distinct"="[r]/[l] contrast")) + # Tweak cosmetics:
  scale_shape(labels=c("yes"="[r] is main allophone", "no"="[r] is not main allophone")) + 
  theme_timo + 
  theme(axis.text.x = element_text(angle = 45, hjust = 1, face = 'bold', size = 10),
        axis.text.y = element_text(face = 'bold', size = 10),
        axis.title = element_text(face = 'bold', size = 12),
        axis.title.y = element_text(face = 'bold', size = 12, margin = margin(t = 0, r = 35, b = 0, l = 0)),
        plot.title = element_text(face = 'bold', size = 14, margin = margin(t = 0, r = 0, b = 30, l = 0)),
        legend.text = element_text(size = 10),
        legend.title = element_blank(),
        legend.position = c(0.95, 0.10),
        legend.justification = c('right', 'bottom'))

# Save:
ggsave(plot = plot_field, filename = "./plots/figure_model_languages_field.pdf", width = 6, height = 6);
ggsave(plot = plot_field, filename = "./plots/figure_model_languages_field.jpg", width = 6, height = 6);
ggsave(plot = plot_field, filename = "./plots/figure_model_languages_field.tif", width = 6, height = 6, compression="lzw", dpi=600);

plot_field;

Figure 45. The distribution of the proportion of matches per language. Please note that this was obtained by fitting a model that includes only the relevant L1 characteristics, Match ~ 1 + r_l_distinction_L1 + trill_real_L1 + trill_occ_L1 + (1 | Language). Figure generated using R version 4.3.3 (2024-02-29)

6.3 Both experiments

p <- gather_draws(b2_res_webfield$b3, `b_.*`, regex = TRUE) %>%
  # remove the b_ for plotting
  #mutate(.variable = gsub("b_","", .variable)) %>%
  ggplot(aes(x = .value, y = .variable, fill = after_stat(x < 0))) +
  geom_vline(xintercept=0.0, linetype="dotted", color="gray30") +
  stat_halfeye(alpha=0.75, point_interval=median_qi, .width=c(.50, .95)) +
  scale_fill_manual(values = c("lightsalmon", "skyblue"), name=NULL, labels=c("≥ 0", "< 0")) +
  xlim(-10, 10) +
  scale_y_discrete(labels=c("b_Intercept" = "Intercept α", 
                            "b_datasetfield" = "the field vs the web dataset",
                            "b_r_l_distinction_L1distinct" = "are [r] and [l] contrasting in L1(s)?",
                            "b_trill_real_L1yes" = "is [r] the main allophone in L1(s)?",
                            "b_r_l_distinction_L2distinct" = "are [r] and [l] contrasting in L2(s)?",
                            "b_trill_real_L2yes" = "is [r] the main allophone in L2(s)?")) +
  ylab(NULL) + xlab("Estimate") + 
  theme_timo + # Tweak cosmetics:
  theme(axis.title.x = element_text(margin = margin(t = 12, b = 0,
                                                    r = 0, l = 0),
                                    face = 'bold', size = 14),
        axis.text.x = element_text(face = 'bold', size = 12),
        axis.text.y = element_text(face = 'bold', size = 12))+
  NULL;
# Show and save:
p;

Figure 46. Posterior estimaes of the fixed effects showing the median (dot), 50% (thick line) and 95% (thin lines) quantiles. Also showing 0 as the vertical dotted line and negative (blue) vs positive (red) values. Figure generated using R version 4.3.3 (2024-02-29)

ggsave(plot = p, filename = "./plots/full_model_both_coefficients.pdf", width = 8, height = 4);
ggsave(plot = p, filename = "./plots/full_model_both_coefficients.jpg", width = 8, height = 4);
ggsave(plot = p, filename = "./plots/full_model_both_coefficients.tif", width = 8, height = 4, compression="lzw", dpi=600);

# take2: generate an exhaustive new dataset:
newdata_webfield$dataset <- factor(c(rep("web", nrow(newdata_web)), rep("field", nrow(newdata_field))), levels=c("web", "field"));

# fit a new model considering only the L1 info:
if( !file.exists("./cached_results/b_plotting_webfield.RData") )
{
  b_plotting_webfield <- brm(Match ~ 1 + dataset + 
                               r_l_distinction_L1 + trill_real_L1 + trill_occ_L1 +
                               (1 | Language) +
                               (1 | Family) +
                               (1 | Autotyp_Area),
                             data=webfield,
                             family=bernoulli(link='logit'),
                             prior=c(brms::set_prior("student_t(5, 0, 2.5)", class="b")),
                             save_pars=save_pars(all=TRUE), # needed for Bayes factors
                             sample_prior=TRUE,  # needed for hypotheses tests
                             seed=998, cores=brms_ncores, iter=10000, warmup=4000, thin=2, control=list(adapt_delta=0.999, max_treedepth=13));
  summary(b_plotting_webfield); mcmc_plot(b_plotting_webfield, type="trace"); mcmc_plot(b_plotting_webfield); # very decent
  bayestestR::hdi(b_plotting_webfield, ci=0.95);
  
  # save the model:
  saveRDS(b_plotting_webfield, "./cached_results/b_plotting_webfield.RData", compress="xz");
} else
{
  b_plotting_webfield <- readRDS("./cached_results/b_plotting_webfield.RData");
}
fit <- fitted(b_plotting_webfield, newdata=newdata_webfield, re_formula=NULL, robust=TRUE);
colnames(fit) = c('fit', 'se', 'lwr', 'upr');
newdata_webfield <- cbind(newdata_webfield, fit);

# Order predictions by descriptive average:
newdata_webfield <- arrange(newdata_webfield, fit);
#newdata_webfield <- mutate(newdata_webfield, Language = factor(Language, levels = newdata_webfield$Language));

# How many languages are over 0.5? (will be reported in paper)
sum(newdata_webfield$lwr > 0.5);

# Finally, add the averages to the plot:
newdata_webfield$avg[ newdata_webfield$dataset == "web" ]   <- web_avg[match(newdata_webfield$Language[ newdata_webfield$dataset == "web" ],     web_avg$Language), ]$M;
newdata_webfield$avg[ newdata_webfield$dataset == "field" ] <- field_avg[match(newdata_webfield$Language[ newdata_webfield$dataset == "field" ], field_avg$Language), ]$M;

# Match language names into there and order:
newdata_webfield$Language[ newdata_webfield$dataset == "web" ]   <- paste0(web[match(newdata_webfield$Language[ newdata_webfield$dataset == "web" ],     web$Language),   ]$Name, " (W)");
newdata_webfield$Language[ newdata_webfield$dataset == "field" ] <- paste0(field[match(newdata_webfield$Language[ newdata_webfield$dataset == "field" ], field$Language), ]$Name, " (F)");
newdata_webfield <- mutate(newdata_webfield, Language = factor(Language, levels = newdata_webfield$Language));

# Setup the plot:
# Aesthetics and geom:
plot_field <- newdata_webfield %>% 
  ggplot(aes(x = Language, col = r_l_distinction_L1, y = fit, ymin = lwr, ymax = upr, shape = trill_real_L1)) +
  #scale_shape_manual(values = c("no" = 1, "yes" = 3)) +
  geom_errorbar(aes(col = r_l_distinction_L1), linewidth = 1.0, width = 0.3) +
  geom_point(size = 5) +
  geom_hline(yintercept = 0.5, linetype = 2, linewidth = 1.5, col = 'grey') + 
  geom_point(aes(y = avg, fill=dataset), col = 'black', shape = 23, size = 5, stroke = 1.0, alpha = 0.25)  +
  ylim(0.5,1.00) +
  labs(x = '', y = 'Proportion of congruent responses') + # Axis labels:
  #ggtitle('Posterior medians and descriptive averages of congruent responses\nby language and R/L contrast (color) and [r] as primary rhotic in L1 (shape)') +
  scale_color_manual(values = c("same"=colorBlindBlack8[2], "distinct"=colorBlindBlack8[3]), labels=c("same"="no [r]/[l] contrast", "distinct"="[r]/[l] contrast")) + # Tweak cosmetics:
  scale_shape(labels=c("yes"="[r] is main allophone", "no"="[r] is not main allophone")) + 
  scale_fill_manual(values = c("web"="white", "field"="black"), labels=c("web"="web experiment", "field"="field experiment")) + 
  theme_timo + 
  theme(axis.text.x = element_text(angle = 45, hjust = 1, face = 'bold', size = 10),
        axis.text.y = element_text(face = 'bold', size = 10),
        axis.title = element_text(face = 'bold', size = 12),
        axis.title.y = element_text(face = 'bold', size = 12, margin = margin(t = 0, r = 35, b = 0, l = 0)),
        plot.title = element_text(face = 'bold', size = 14, margin = margin(t = 0, r = 0, b = 30, l = 0)),
        legend.text = element_text(size = 10),
        legend.title = element_blank(),
        legend.position = c(0.95, 0.10),
        legend.justification = c('right', 'bottom'))

# Save:
ggsave(plot = plot_field, filename = "./plots/figure_model_languages_field.pdf", width = 18, height = 6);
ggsave(plot = plot_field, filename = "./plots/figure_model_languages_field.jpg", width = 18, height = 6);
ggsave(plot = plot_field, filename = "./plots/figure_model_languages_field.tif", width = 18, height = 6, compression="lzw", dpi=600);

plot_field;

Figure 47. The distribution of the proportion of matches per language. Please note that this was obtained by fitting a model that includes only the relevant L1 characteristics, Match ~ 1 + dataset + r_l_distinction_L1 + trill_real_L1 + trill_occ_L1 + (1 | Language) + (1 | Family) + (1 | Autotyp_Area). Figure generated using R version 4.3.3 (2024-02-29)

7 Conclusions

It is clear that across the two datasets our participants have a strong tendency to associate [r] with the zigzagy line and [l] with the straight line well above the conservative 50% chance level. Presenting [r] first seems to increase the probability of a match in the online experiment, compared to presenting [l] first. The field experiment has an overall higher probability of a match than the online experiment. In the online experiment and in the combined datasets, having [r] as the primary allophone in the the L1(s) spoken by a participant slightly decreases the probability of a match.

8 Session information

if( require(benchmarkme) )
{
  # CPU:
  cpu_info <- benchmarkme::get_cpu();
  if( is.null(cpu_info) || all(is.na(cpu_info)) )
  {
    cat("**CPU:** unknown.\n\n");
  } else
  {
    if( !is.null(cpu_info$model_name) && !is.na(cpu_info$model_name) )
    {
      cat(paste0("**CPU:** ",cpu_info$model_name));
      if( !is.null(cpu_info$no_of_cores) && !is.na(cpu_info$no_of_cores) )
      {
        cat(paste0(" (",cpu_info$no_of_cores," threads)"));
      }
      cat("\n\n");
    } else
    {
      cat("**CPU:** unknown.\n\n");
    }
  }
  
  # RAM:
  ram_info <- benchmarkme::get_ram();
  if( is.null(ram_info) || is.na(ram_info) )
  {
    cat("**RAM (memory):** unknown.\n\n");
  } else
  {
    cat("**RAM (memory):** "); print(ram_info); cat("\n");
  }
} else
{
  cat("**RAM (memory):** cannot get info (try installing package 'benchmarkme').\n\n");
}

CPU: Apple M3 (8 threads)

RAM (memory): 17.2 GB

pander::pander(sessionInfo());

R version 4.3.3 (2024-02-29)

Platform: aarch64-apple-darwin20 (64-bit)

locale: en_US.UTF-8||en_US.UTF-8||en_US.UTF-8||C||en_US.UTF-8||en_US.UTF-8

attached base packages: parallel, stats, graphics, grDevices, utils, datasets, methods and base

other attached packages: benchmarkme(v.1.0.8), knitr(v.1.48), sjPlot(v.2.8.16), ggpubr(v.0.6.0), gridExtra(v.2.3), bayestestR(v.0.14.0), performance(v.0.12.2), lme4(v.1.1-35.5), Matrix(v.1.6-5), tidybayes(v.3.0.6), ggdist(v.3.3.2), brms(v.2.21.0), Rcpp(v.1.0.13), lubridate(v.1.9.3), forcats(v.1.0.0), stringr(v.1.5.1), dplyr(v.1.1.4), purrr(v.1.0.2), readr(v.2.1.5), tidyr(v.1.3.1), tibble(v.3.2.1), ggplot2(v.3.5.1) and tidyverse(v.2.0.0)

loaded via a namespace (and not attached): RColorBrewer(v.1.1-3), tensorA(v.0.36.2.1), rstudioapi(v.0.16.0), jsonlite(v.1.8.8), datawizard(v.0.10.0), magrittr(v.2.0.3), TH.data(v.1.1-2), estimability(v.1.5.1), farver(v.2.1.2), nloptr(v.2.1.1), rmarkdown(v.2.27), ragg(v.1.3.2), vctrs(v.0.6.5), minqa(v.1.2.7), rstatix(v.0.7.2), htmltools(v.0.5.8.1), distributional(v.0.4.0), curl(v.5.2.1), haven(v.2.5.4), broom(v.1.0.6), sjmisc(v.2.8.10), sass(v.0.4.9), StanHeaders(v.2.32.9), bslib(v.0.8.0), plyr(v.1.8.9), sandwich(v.3.1-0), emmeans(v.1.10.2), zoo(v.1.8-12), cachem(v.1.1.0), TMB(v.1.9.14), iterators(v.1.0.14), lifecycle(v.1.0.4), pkgconfig(v.2.0.3), sjlabelled(v.1.2.0), R6(v.2.5.1), fastmap(v.1.2.0), digest(v.0.6.36), numDeriv(v.2016.8-1.1), colorspace(v.2.1-1), textshaping(v.0.4.0), labeling(v.0.4.3), fansi(v.1.0.6), timechange(v.0.3.0), httr(v.1.4.7), abind(v.1.4-5), mgcv(v.1.9-1), compiler(v.4.3.3), pander(v.0.6.5), doParallel(v.1.0.17), bit64(v.4.0.5), withr(v.3.0.1), backports(v.1.5.0), inline(v.0.3.19), carData(v.3.0-5), QuickJSR(v.1.3.1), pkgbuild(v.1.4.4), highr(v.0.11), ggsignif(v.0.6.4), MASS(v.7.3-60.0.1), sjstats(v.0.19.0), loo(v.2.7.0), tools(v.4.3.3), glue(v.1.7.0), nlme(v.3.1-165), grid(v.4.3.3), checkmate(v.2.3.2), reshape2(v.1.4.4), generics(v.0.1.3), glmmTMB(v.1.1.9), gtable(v.0.3.5), tzdb(v.0.4.0), hms(v.1.1.3), car(v.3.1-2), utf8(v.1.2.4), foreach(v.1.5.2), pillar(v.1.9.0), vroom(v.1.6.5), posterior(v.1.5.0), benchmarkmeData(v.1.0.4), splines(v.4.3.3), lattice(v.0.22-6), survival(v.3.7-0), bit(v.4.0.5), tidyselect(v.1.2.1), arrayhelpers(v.1.1-0), V8(v.4.4.2), stats4(v.4.3.3), xfun(v.0.46), bridgesampling(v.1.1-2), matrixStats(v.1.3.0), rstan(v.2.32.6), stringi(v.1.8.4), yaml(v.2.3.10), boot(v.1.3-30), evaluate(v.0.24.0), codetools(v.0.2-20), cli(v.3.6.3), RcppParallel(v.5.1.8), systemfonts(v.1.1.0), xtable(v.1.8-4), munsell(v.0.5.1), jquerylib(v.0.1.4), ggeffects(v.1.7.0), coda(v.0.19-4.1), svUnit(v.1.0.6), rstantools(v.2.4.0), bayesplot(v.1.11.1), Brobdingnag(v.1.2-9), viridisLite(v.0.4.2), mvtnorm(v.1.2-5), scales(v.1.3.0), insight(v.0.20.2), crayon(v.1.5.3), rlang(v.1.1.4) and multcomp(v.1.4-26)